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Introduction 


Background 

This  study  relates  to  the  further  development  of  MILHY.  Two 
previous  reports  are  relevant  to  the  research  reported  here. 
In  the  first  of  these  two  reports  (1),  a  review  of  available 
hydrological  models  was  undertaken,  and  a  case  was  made  for 
the  further  development  of  MILHY  as  an  operational  model  for 
ungauged  catchment  flood  forecasting.  The  subsequent  report 
(2)  detailed  two  applications  of  a  revised  MILHY  scheme 
(referred  to  here  as  MILHY2),  in  which  the  curve  number 
scheme  for  the  estimation  of  runoff  was  replaced  by  a  finite 
difference  scheme.  The  advantage  of  such  a  replacement  was 
seen  to  be  the  improved  time  resolution  of  runoff 
prediction  and  the  improved  accommodation  of  anticedent 
conditions  whilst  retaining  the  same  data  input  requirements 
as  MILHY.  The  results  of  the  two  applications  undertaken 
were  sufficiently  encouraging  for  the  model  development  work 
to  be  continued,  and  it  is  this  work  that  is  the  subject  of 
the  current  report. 

Obiectlves  and  Scone 


The  two  principal  objectives  of  the  research  work  reported 
here  were: 


(i)  The  application  of  MILHY2  to  further  watersheds, 
(ii)  The  presentation  of  the  Fortran  program  for  MILHY2 


Figure  1  illustrates  how  these  objectives  fit  into  the 


author's  view  of  the  conceptual  and  operational  development 
of  MILHY2,  as  outlined  at  the  MILHY  Workshop  st  W.E.S.  on 
the  12  January  1985.  In  that  outline,  the  work  reported 
here,  and  the  objectives  above,  relate  to  the  Increase  in 
validation  (operational)  and  to  the  development  of  the 
Fortran  version  of  MILHY2  (conceptual)  under  1985. 
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Figure  1  :  Conceptual  and  operational  ^'^'-elopments 
of  MILHY  research 


Further  Application  of  MILHY2 


2 . I  Introduction 

Certain  results  of  the  application  of  MILHY2  to  the  North  Creek 
catchment,  Texas  and  the  Sixmile  Creek  catchment,  Arkansas,  have  been 
presented  in  DAJA37-81-C-0221  in  the  context  of  operational  validation. 
Application  to  these  catchments  was  used  firstly  to  Illustrate  the 
suitability  of  the  Brakensiek  and  Rawls  empirical  information  for  the 
derivation  of  the  soils  data  necessary  to  operate  the  model  quite 
successfully  for  the  ungauged  catchment,  secondly  to  illustrate  a 
favourable  comparison  of  calculated  to  measured  hydrographs  for  certain 
experimental  frames.  The  deterministic  version  of  MILHY2  is  thus 
considered  to  be  operationally  valid  for  the  variety  of  conditions  which 
have  been  considered  so  far.  .However,  it  is  important  to  extend  this 
range  of  application  and  consequently,  the  details  of  the  application  of 
MILHY2  to  a  further  five  catchments  in  Vermont  and  Iowa,  United  States 
of  America  are  provided  in  this  report  together  with  program  code.  In 
addition,  these  applications  will  provide  information  for  discussion  of 
the  following  points: 

1  Is  MILHY2  of  a  form  which  is  suitable  for  application  to  the  ungauged 
catchment  ? 

The  runoff  procedure  which  has  been  introduced  in  DAJA37-8 l-C-022 1  is 
not  a  simple  calibrated  procedure,  but  is  physically  based.  Much  of 
the  original,  and  so  far  undeveloped,  model  however,  does  remain 
calibrated  and  the  issue  of  the  validity  of  extrapolation  of  results 
which  have  been  produced  by  calibration  to  other  gauged  catchments 


must  be  raised. 

2  Can  MILHY2  meet  an  operational  requirement? 

Operational  requirements  were  discussed  in  (2).  It  has  already  been 
established  that  MILHY2  can  be  ported  onto  a  microcomputer  system. 
Application  will  reveal  whether  or  not  the  model  will  run  at 
acceptable  speeds  on  this  hardware  configuration.  In  addition,  the 
following  questions  must  be  considered: 

a  Are  the  data  preparation  requirements  reasonable  in  the  context  of  a 
potential  nonprofessional  user? 

b  Can  sufficient  guidelines  be  provided  for  the  user  in  terms  of 
application  and  interpretation  of  the  model  for  a  range  of 
applications  ? 

c  Can  the  model  be  made  user  friendly? 

d  Is  the  software  reliable  for  the  now  expanding  range  of  applications? 

3  Does  MILHY2  have  an  appropriate  structure  for  the  ungauged  and 
operational  application? 

The  physically  based  infiltration  model  which  has  been  developed, 
although  simple,  does  attempt  to  attain  a  balance  between  a 
methodology  which  is  scientifically  acceptable,  and  one  which  remains 
operationally  feasible.  The  suitability  of  this  choice  will  be 
revealed  with  the  application  of  MILHY2. 

In  any  application,  there  will  be  interest  in  the  accuracy  of  the 
hydrograph  predictions  which  the  model  supplies.  However,  it  has  been 
stressed  throughout  the  discussion  on  model  evaluation,  that  there  are 
other  important  questions  which  must  also  be  specifically  investigated 
in  order  to  provide  an  unskilled  user  with  sufficient  information  to 
guide  the  intelligent  use  of  the  model.  In  addition  to  a  comparison  of 
calculated  and  measured  hydrographs,  the  following  questions  must  also 
be  addressed  during  application  of  MILHY2: 

1  What  is  involved  in  the  data  acquisition  and  preparation  stage?  A 
user  needs  to  know  the  nature  of  the  decisions  which  must  be  taken  in 
order  to  derive  the  necessary  model  parameters.  It  is  also  important 


to  assess  the  likely  time  period  which  will  be  required  for  data 
preparation. 

2  Is  the  infiltration  behaviour  predicted  by  the  physically  based 
infiltration  model  reasonable  for  a  range  of  catchment  situations? 
Infiltration  behaviour  has  been  examined  for  a  range  of  hypothetical 
conditions.  It  is  important  to  examine  its  behaviour  for  more 
complex  soil  and  precipitation  conditions. 

3  Is  the  explicit  finite  difference  method  accurate  for  these  more 
complex  soil  profile  and  variable  storm  conditions? 

These  issues  are  now  considered  specifically  for  catchment  situations. 
These  three  issues:  data  preparation,  infiltration  behaviour,  and  the 
stability  of  the  numerical  solution,  have  not  been  discussed  in  the 
context  of  the  application  to  the  North  Creek  and  Sixmile  Creek 
catchments.  The  information  derived  from  these  two  catchments  will 
therefore  be  included  in  those  relevant  sections. 

This  report  will  therefore  be  divided  into  six  sections.  Firstly,  the 
five  catchments  which  are  to  be  used  in  this  chapter  will  be  introduced 
(2.2).  Secondly,  the  data  collection  and  preparation  which  are 
necessary  for  the  application  of  MILHY2  to  the  catchments  will  be 
described  (2.3).  In  addition,  some  more  general  points  about  this 
critical  stage  in  model  application  will  be  made.  Thirdly,  a  series  of 
comparisons  of  calculated  and  measured  hydrographs  for  a  range  of 
storms,  applied  to  the  five  catchments  in  Vermont  and  Iowa,  will  be 
presented  (2.4).  This  comparison  will  follow  the  two  stage  procedure  in 
figure  2.  Fourthly,  the  infiltration  behaviour  which  is  predicted  by 
the  model  for  the  layered  soil  profiles  and  more  erratic  rainfall 
conditions,  experienced  by  the  catchments  and  the  numerical  errors 
incurred  in  the  solution  of  the  Richards  equation  by  the  explicit  finite 
difference  method  will  be  examined  (section  3).  Finally,  an  attempt 
will  be  made  to  summarize  the  information  derived  from  all  experimental 
frames  which  have  been  used,  in  order  to  define  those  conditions  for 
which  the  model  is,  and  those  for  which  it  is  not,  appropriate  (section 
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2.2  Catchment  location  details 


The  five  catchments  documented  In  this  chapter,  and  which  have  been  used 
to  evaluate  the  operation  of  MILHY2  are  the  following: 

1  An  unnamed  tributary  of  the  Sleepers  River  catchment,  Connecticut 
River  basin,  watershed  2  (W-2)  in  North  Danville,  Vermont,  United 
States  of  America. 

2  Watershed  I  (W-l),  Silver  Creek,  West  Nlshnabotna  River,  Missouri 
River  basin,  Treynor,  Iowa,  United  States  of  America. 

3  Watershed  2  (W-2),  Keg  Creek,  Missouri  River  basin,  Treynor,  Iowa, 
United  States  of  America. 

4  Watershed  3  (W-3),  Silver  Creek,  West  Nlshnabotna  River,  Missouri 
River  basin,  Treynor,  Iowa,  United  States  of  America. 

3  Watershed  4  (W-4),  Silver  Creek,  West  Nlshnabotna  River,  Missouri 
River  ba.iln,  Treynor,  Iowa,  United  States  of  America, 

The  Location  of  these  catchments  is  indicated  in  figure  3,  and  a 
comparison  of  the  three  physical  catchment  characteristics  which  are 
required  by  the  unit  hydrograph  procedure.  Is  provided  by  table  1.  All 
of  these  catchments  are  small  in  area  (less  than  0.6  square  km)  as  this 
enables  a  closer  examination  to  be  made  of  the  modified  runoff  component 
of  the  model  wit!iout  incorporating  the  need  for  channel  routing. 

All  of  these  catchments  are  gauged  catchments  and  are  United  States 
Department  of  Agriculture  (USDA)  Agricultural  Research  Service  (ARS) 
experimental  watersheds.  Hydrological  data  from  all  ARS  experimental 
watersheds  ire  currently  stored  on  a  data  base  in  the  United  States, 
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0.6 
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which  is  accessible  by  use  of  REPHLEX  (REtrieval  Procedures  for 
HydroLogic  data  from  ARS  Experimental  watersheds)  which  has  been 
developed  by  the  Water  Data  Laboratory  and  documented  by  Thurman  et  al 
(3).  This  data  base  provides  information  for  305  watersheds  which  range 
from  0.2  ha  to  536  square  km  in  area.  Precipitation  and  runoff  data  for 
individual  storm  events  and  for  daily,  monthly,  or  annual  accumulations, 
and  which  range  in  length  of  record  from  1  to  45  years  are  available. 
Information  may  be  derived  from  the  system  in  tabular  or  graphical  form. 
An  inventory  of  the  ARS  experimental  watersheds  (4)  is  published  which 
documents  the  types  of  data  (precipitation,  runoff,  pan  evaporation, 
soil  moisture,  land  use,  soil  survey,  for  example)  which  are  available 
for  each  catchment. 

The  Sleepers  River  catchment,  Connecticut  River  basin,  Vermont,  is 
located  8.05  km  north  west  of  St.  Johnsbury.  This  catchment  has  been 
the  location  of  many  field  studies  including  Dunne  and  Black  (5,6)  and 
it  is  considered  to  represent  a  typical  glaciated  upland  catchment  of 
New  England.  The  location  and  physical  characteristics  of  the  unnamed 
tributary  W-2,  are  indicated  in  figure  4.  It  is  described  by  the  USDA 
as  comprising  sloping  to  steep  land  at  higher  elevations.  It  has  a 
covering  of  glacial  till  which  exhibits  good  surface  drainage  and  which 
overlies  Devonian  schist  interbedded  with  limestone.  The  land  use 
within  the  watershed  W-2  is  divided  between  permanent  hay  (37X),  pasture 
(38%),  and  maple  and  beech  trees  (25%). 

The  four  catchments  near  Treynor,  Iowa  contain  soils  which  have 
developed  from  the  deep  mantle  of  Wisconsin  loess  (3.05  to  27.72  metres) 
which  overlies  Kansan  glacial  till  which  in  turn  overlies  the  bedrock  of 
interbedded  calcareous  shales  and  limestones.  The  watershed  topography 
has  developed  totally  by  erosion  of  loess  and  the  deeper  gullies  have 
incised  slightly  into  the  till.  The  loess  is  considered  to  have  a 
moderate  rate  of  percolation.  In  all  four  watersheds  channel  flow  is 
permanent  and  fed  by  a  zone  of  saturation  and  seepage  which  occurs  at 
the  loess  and  till  Interface.  Topographic  maps  of  the  four  catchments 
are  provided  in  figures  5-8.  W-1  is  located  9.65  km  south  west  of 
Treynor.  The  catchment  is  laid  to  contour  corn  and  exhibits  high  levels 
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Watershed  3,  Silver  Creek,  West  Nishnabotna  River,  Missouri 
River  Basin,  Treynor,  Iowa 


of  fertility  and  good  farming  practices.  W-2,  also  9.65  km  south  west 
of  Trynor,  has  similar  characteristics  to  W-1  but  is  a  tributary  of 
another  stream,  the  Keg  Creek.  W-3  is  located  4.83  km  south  west  of 
Treynor  and  contains  pasture  with  controlled  grazing.  Finally,  W-4, 
located  4.83  km  south  west  of  Treynor,  contains  contour  corn  on  grassed 
backed  slope  terraces.  All  terraces  in  W-4  are  as  recommended  by  the 
ARS. 

The  five  catchments  which  have  been  introduced  here  are  all  below  0.6 
square  km.  Although  these  may  be  considered  to  be  small,  certain 
limitations  are  imposed  upon  the  catchment  scale  by  the  nature  of  a 
three  year  research  programme.  Within  a  three  year  period,  it  is 
considered  that  three  potential  research  strategies  are  feasible  within 
a  geographical  hydrological  modelling  exercise. 

Firstly,  at  one  extreme,  it  would  be  possible  to  develop  and  implement 
an  entirely  new  mathematical  hydrological  model.  This  would  demand  such 
an  Investment  of  time  that  evaluation  and  testing  could  only  be 
undertaken  for  one  catchment.  Secondly,  it  would  be  possible  to  provide 
a  modification  to  one  component  of  a  currently  utilized  hydrological 
model,  thus  allowing  sufficient  time  for  a  more  detailed  evaluation  of 
the  modified  model  on  a  series  of  catchments  exhibiting  different 
characteristics.  Thirdly,  and  at  the  other  extreme,  it  would  be 
possible  to  apply  a  currently  used  model  to  a  very  large  number  of 
catchments,  but  to  provide  no  model  development.  In  this  third 
strategy,  a  broader  and  more  comprehensive  model  evaluation  could  be 
accomplished . 

The  first  strategy  has  been  a  very  popular  choice.  Fenves  et  al  (7) 
stressed  that  emphasis  has  been  placed  upon  model  development  whilst 
support,  documentation,  and  evaluation  have  been  neglected.  This  has 
led  to  a  multiplicity  of  mostly  underutilized  models  with  no  clear 
recommendations  for  future  requirements  and  research.  Certainly  durirfg 
a  three  year  research  period.  Insufficient  time  would  remain  after  model 
design  and  implementation  fully  to  evaluate  the  model  and  to  examine  its 
full  potential. 


The  third  strategy  has,  in  comparison,  not  commonly  been  undertaken.  It 
has  been  stressed  that  model  evaluation  has  not  been  a  popular 
occupation  in  mathematical  hydrological  modelling.  However,  although 
providing  an  opportunity  for  a  comprehensive  model  evaluation  and 
examination  of  operational  applications,  the  third  strategy  would  not 
allow  for  an  investigation  of  ungauged  catchment  applications  as  no 
suitable  model  could  be  identified.  It  would  also  not  allow  for  the 
examination  of  the  potential  of  a  physically  based,  rather  than  an 
empirical  model  for  application  purposes. 

These  issues  were  considered  to  be  of  importance  and  therefore  the 
second  strategy  was  adopted  in  this  analysis.  A  modification  to  the 
Infiltration  component  of  HYMO  was  undertaken,  and  the  period  of  model 
modification  and  Implementation  has  necessarily  limited  the  available 
time  for  catchment  selection,  data  collection,  and  preparation.  Thus 
seven  small  catchments  were  chosen.  This  provides  a  good  compromise 
between  the  time  limitations  of  a  three  year  research  programme  and  the 
need  to  evaluate  the  model  over  a  range  of  catchment  conditions. 

The  small  size  of  catchments  is  not  a  disadvantage  because  the  emphasis 
in  this  investigation  of  HYMO  and  MILHY2  has  concentrated  upon  the 
hydrograph  computation  procedure.  It  has  not  been  designed  to  examine 
the  characteristics  of  the  Variable  Storage  Coefficient  channel  routing 
technique.  The  selection  of  smaller  catchments  which  can  in  the  context 
of  the  application  of  MILHY2  be  treated  as  single  subcatchments,  has 
allowed  the  hydrograph  computation  to  be  investigated  without  the 
complications  of  the  incorporation  of  the  routing  procedure. 
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The  five  catchments  which  have  now  been  Introduced  in  section  2.2  are 
all  less  than  0.6  square  km  (table  1).  No  subdivision  of  catchments  has 
been  necessary,  and  consequently  no  channel  cross  section  information  is 
required  for  channel  routing  operations.  The  catchment  characteristics: 
area,  elevation  difference,  and  main  stream  length  (table  1),  have  been 
derived  for  all  five  catchments  from  maps  of  the  scale  and  detail 
illustrated  in  figures  4  to  8.  The  determination  of  the  soils  data  will 
now  be  discussed  for  each  catchment. 

There  are  five  major  soil  types  in  the  Watershed  W-2,  North  Danville, 
Vermont.  These  Include  sandy  loams,  silt  loams,  and  loams,  and  are 
namely,  Colraln,  Peacham,  Calais,  Cabot,  and  Woodstock.  The  details 
concerning  soil  horizon  depths  and  soil  textural  characteristics  of  each 
layer  were  available  from  the  USDA  ARS  descriptions  of  the  catchment 
(table  2).  The  division  of  each  soil  horizon  into  cells  was 
accomplished  according  to  the  general  rule  that  cells  in  the  top  layer 
must  not  be  greater  than  0.1  metres  and  in  the  lower  two  layers,  not 
greater  than  0.15  metres.  From  the  soil  texture  information,  the 
Brakensiek  and  Rawls  charts  were  used  to  define  the  soil  hydrological 
characteristics.  For  all  soil  textures,  the  centroid  position  on  the 
Bakensiek  and  Rawls  charts  was  used.  Detention  capacity  was  assumed  to 
be  zero  and  a  uniform  initial  relative  saturation  of  80%  was  assumed. 

The  four  catchments  near  Treynor  all  contain  the  same  four  soil  types, 
but  each  soil  occupies  different  proportions  of  the  total  catchment 
area.  The  four  soil  types  are  Monona,  Marshall,  Napier,  and  Ida,  and 
comprise  silt  loams  and  silty  clay  loams.  Very  little  information  was 
available  on  the  layering  characteristics  of  these  soils  and  therefore, 
no  layering  of  the  representive  soil  columns  was  Incorporated.  The 
hydrological  characteristics  of  each  soil  texture  group  were  derived 
from  the  centroid  position  of  the  Brakensiek  and  Rawls  charts.  The  soil 
column  which  is  defined  by  the  depth  of  the  soil  is  divided  into  equal 
sized  cells  of  0.05  metres  for  Napier  (the  deepest  soil)  and  0.025 


metres  for  the  other  three  soils.  The  details  of  the  soils  in  all  four 
of  these  catchments  are  provided  by  table  2.  The  detention  capacity  of 
catchment  W-4  was  set  at  0.01  metres.  This  value  is  estimated  according 
to  the  terracing.  No  detention  capacity  was  assumed  for  the  other  three 
catchments.  Initial  relative  saturation  was,  in  the  absence  of  soil 
moisture  information  and  based  on  previous  experience,  assumed  to  be  80% 
at  the  surface,  and  to  Increase  uniformly  with  depth. 

The  precipitation  data  for  all  storms  applied  to  these  five  catchments 
were  converted  into  cumulative  totals  at  equal  time  intervals,  the  form 
which  is  required  by  MILHY2.  The  measured  hydrograph  for  each  storm 
event  was  also  input  to  MILHY2  for  comparison.  The  storms  which  were 
used  and  the  runoff  which  they  produced  are  indicated  in  table  3. 

Experience  of  application  of  the  model  to  these  five  catchments,  and 
those  of  Texas  and  Arkansas,  has  illustrated  that  in  order  to  provide 
the  data  for  model  application,  the  user  is  involved  in  four  stages. 
Figure  9  illustrates  these  stages,  which  include  data  collection,  data 
preparation,  data  entry  and  data  checking. 

Data  collection 

This  involves  securing  three  sources  of  information:  a  topography  map  of 
the  catchment,  a  soils  map  and  accompanying  description,  and 
precipitation  data.  Depending  upon  the  level  of  information  which  is 
available,  the  precipitation  data  might  be  in  the  form  of  recording  rain 
gauge  data,  storm  totals  or  predicted  rainfall  data.  The  distribution 
of  precipitation,  where  only  storm  totals  are  available,  may  be  provided 
by  application  of  one  of  the  standard  Soil  Conservation  Service  rainfall 
distribution  models. 

Data  preparation 

This  involves  the  user  in  a  number  of  decisions  as  to  the  manner  in 
which  the  catchment  should  be  characterized,  the  use  of  the  Brakensiek 
and  Rawls  tables  and  cnarts  to  derive  soil  hydrological  properties,  and 
a  series  of  manual  calculations  to  convert  precipitation  data  into  the 
form  required  by  M1LHY2.  All  of  these  actions  could  potentially 


Soil  type  USDA  Average  Catchment 

soil  texture  depth  of  area 
soil  (%) 

(metres) 


W-2  North  Danville 

,  Vermont 

Colrain 

sandy  loam 

0.84 

41 

Peacham 

silt  loam 

0.31 

5 

Calais 

loam 

0.69 

9 

Cabot 

silt  loam 

0.46 

13 

Woodstock 

sandy  loam 

0.61 

32 

Treynor ,  Iowa 

W-1 

W-2 

W-3 

W-4 

Monona 

silt  loam 

0.15 

38 

24 

50 

48 

Marshall 

silty  clay  loam 

0.254 

35 

36 

22 

23 

Napier 

silt  loam 

0.762 

16 

17 

22 

23 

Ida 

silt  loam 

0.076 

11 

23 

6 

6 

Storm 

Date  of 

Time  of 

Time 

Storm 

Total 

Total 

number 

storm  start 

storm 

increment 

duration 

precip¬ 

runoff 

(d.m.yr) 

start 

of  rainfall 

(hrs) 

itation 

(mm) 

(hrs)  (hrs)  (mm) 


W-2, 

North  Danville, 

Vermont 

1 

11.9.1968 

06:00 

1.0 

16.0 

38.1 

0.36 

2 

21.7.1969 

15.30 

0.25 

3.25 

24.1 

0.31 

3 

28.8.1970 

14:45 

0.25 

6.5 

37.3 

0.54 

4 

16.7.1967 

04:30 

0.5 

9.0 

43.9 

4.67 

5 

30.7.1960 

12:00 

1.0 

11.0 

43.9 

2.72 

6 

2.6.1961 

02:00 

0.25 

6.0 

21.1 

4.39 

W-1, 

Treynor,  Iowa 

1 

2.8.1970 

21:40 

0.1 

1.8 

67.1 

22.96 

2 

26.6.1966 

02:32 

0.1 

1.0 

22.9 

9.27 

3 

14.6.1967 

05:10 

0.1 

1.7 

19.6 

12.34 

4 

20.6.1967 

20:56 

0.05 

2.9 

156.0 

107.30 

5 

7.6.1967 

17:05 

0.1 

1.4 

41.9 

31.3 

W-2, 

Treynor,  Iowa 

1 

2.8.1970 

21:37 

0.1 

1.8 

41.9 

17.96 

2 

26.6.1966 

02:26 

*0.1 

1.2 

22.9 

10.19 

3 

14.6.1967 

05:13 

0.1 

1.7 

19.8 

10.97 

4 

20.6.1967 

20:56 

0.05 

2.75 

143.0 

96.16 

5 

7.6.1967 

17:10 

0.1 

1.0 

43.2 

25.62 

W-3, 

Treynor,  Iowa 

1 

2.8.1970 

21:33 

0.1 

1.7 

41.7 

1.52 

2 

25.6.1966 

23:05 

0.1 

1.3 

28.7 

4.l'4 

3 

14.6.1967 

05:10 

0.1 

1.8 

21.1 

2.99 

4 

20.6.1967 

20:52 

0.1 

2.8 

98.6 

33.75 

5 

7.6.1967 

17:10 

0.1 

1.3 

23.9 

4.17 

W-4, 

Treynor,  Iowa 

1 

2.8.1970 

21:33 

0.1 

1.7 

41.7 

0.15 

2 

26.6.1966 

23:05 

0.1 

1.3 

28.7 

1.27 

3 

14.6.1967 

05:10 

o.a 

1.8 

21.1 

1.21 

4 

20.6.1967 

20:52 

0.1 

2.8 

98.6 

9.53 

5 

7.6.1967 

17:10 

0.1 

1.3 

23.9 

1.44 

Data  collection 


Figure  9; 


The  four  stages  in  data  generation 


Introduce  error  into  the  predictions.  To  reduce  this  source  of  error, 
and  to  operationalize  the  model  as  fully  as  possible  for  the 
nonprofessional  hydrologist,  it  is  important  that  an  attempt  should  be 
made  to  computerize  certain  procedures  in  this  data  preparation  stage. 

It  is  necessary  that  the  catchment  characteristics  required  by  the  unit 
hydrograph  method:  catchment  area,  main  stream  length,  and  difference  in 
elevation,  he  determined  by  the  user.  This  is  a  straightforward,  but 
tedious  procedure,  which  does  not  require  specialized  skills.  The 
determination  of  area  could  only  be  computerized  should  a  digitizing 
facility  be  available  on  the  computer  system.  Access  to  this  cannot  be 
assumed  for  the  microcomputer  system  user.  However,  it  is  important  to 
stress  to  the  user  the  Importance  of  accuracy  in  the  specification  of 
these  three  catchment  characteristics.  Figure  10  provides  a  summary  of 
certain  results  of  the  application  of  a  deterministic  sensitivity 
analysis  to  the  unit  hydrograph  method  which  is  used  by  MILHY2.  The 
sensitivity  of  the  peak  unit  discharge  to  the  three  catchment 
characteristics  is  illustrated.  For  a  constant  elevation  difference  of 
15.24  metres,  figure  10(A)  illustrates  that  as  the  area  of  the  catchment 
increases,  i.e.  topography  becomes  less  steep,  the  sensitivity  of  unit 
peak  to  length  of  main  channel  increases.  For  any  given  area  and  height 
combination,  the  sensitivity  to  length  of  main  channel  is  greatest  when 
the  channel  is  shorter.  Figure  10(B)  illustrates  that  the  unit  peak  is 
sensitive  to  catchment  area.  This  sensitivity  is  greatest  for  smaller 
catchment  areas  and  varies  quite  significantly  according  to  the  height 
to  length  ratio.  As  this  ratio  decreases  and  topography  becomes  less 
steep,  then* sensitivity  to  area  increases.  Figure  10(C)  illustrates 
that  the  sensitivity  of  the  model  to  elevation  difference  decreases  as 
the  height  difference  Increases,  The  magnitude  of  this  sensitivity  is 
related  to  the  catchment  shape,  being  less  for  narrower  and  elongated 
catchments.  It  is  important  therefore,  that  these  three  catchment 
characteristics  are  specified  as  accurately  as  possible. 

The  selection  of  the  major  soil  types  is  another  choice  for  which  very 
little  direct  help  can  be  provided  specifically  for  the  catchment  of 
interest  to  the  user.  Examination  of  the  soils  map  is  necessary  to 
Identify  the  major  soil  types,  and  to  determine  the  percentage  of  the 
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catchment  area  which  each  occupies. 

It  is  intended  that  the  experience  of  a  series  of  applications  of 
MILHY2,  which  are  documented  in  this  report,  will  be  useful  in  defining 
a  very  general  series  of  guidelines  to  which  the  user  may  refer  when 
selecting  the  appropriate  number  of  soil  columns  to  represent  the 
catchment  area,  the  layering  characteristics  of  each  soil  column,  and 
the  dimensions  of  the  cells  in  the  soil  column. 

The  number  of  soil  columns  will  reflect  a  trade-off  between  a  possible 
Increase  in  prediction  accuracy  and  the  increased  computer  and  data 
preparation  costs  which  are  associated  with  the  application  of  a  large 
number  of  soil  columns.  If  sufficient  detail  is  available  in  the  soil 
map  descriptions  to  define  Che  soil  texture  characteristic  of  up  to 
three  layers  in  the  soil,  then  this  information  can  be  used.  Should 
this  degree  of  data  not  be  available,  the  user  must  have  access  to 
advice  or  a  standard  procedure  which  can  be  applied.  Choice  of  the  size 

and  hence  Che  number  of  cells  in  the  soil  column  should  also  be  based  on 

the  past  experience  of  application  of  the  model. 

If  a  general  series  of  rules  based  upon  the  results  of  gauged 
applications  on  the  model  caji  indeed  be  established,  then  it  is 
Important  that  a  user  does  have  access  to  this  information.  There  are 
two  forms  in  which  this  information  may  be  stored.  Firstly,  it  can  be 
provided  in  a  manual  which  accompanies  the  computer  program,  or 
secondly,  it  can  be  provided  on-line.  The  Information  can  be  held  in 
the  computer  program  and  provided  to  the  user  on  request,  in  an 

interactive  form,  as  the  user  enters  the  data  for  model  application. 

For  example,  where  the  user  is  required  to  specify  the  number  of  soil 
columns  for  the  catchment  area,  if  insufficient  Information  is 
available,  or  if  the  user  is  unfamiliar  with  the  model,  then  the  user 
may  interrogate  the  system  for  advice.  Based  on  past  application,  Che 
number  of  soil  columns  can  be  related  to  catchment  size,  precipitation 
characteristics,  the  size  of  the  computer  system,  and  to  any  constraints 
which  Che  user  might  be  placing  on  response  time.  The  user  will  then  be 
in  a  position  to  operate  the  model  to  a  greater  advantage  and  based  upon 


the  past  experience  of  the  model  application,  rather  than  on  past 
personal  experience.  With  time,  the  information  which  is  held  by  the 
system  can  be  increased. 

The  use  of  the  Brakensiek  and  Rawls  (see  pages  32,  33  in  (2))  charts  to 
provide  the  soil  hydrological  characteristics,  saturated  hydraulic 
conductivity,  saturated  moisture  content,  and  soil  moisture 
characteristic  curve,  is  one  very  obvious  area  where  operator  error  may 
be  reduced.  The  look-up  procedure  which  uses  the  tables  could  be 
replaced  by  a  series  of  expressions  which  are  more  easily  computerized. 
It  is  only  necessary  for  the  user  to  define  the  soil  texture  class,  sand 
or  loam  for  example,  for  each  soil  type,  and  each  layer  where 
appropriate.  This  information  is  then  entered  into  a  program  which  will 
firstly  convert  the  soil  texture  category  to  a  percentage  clay  and 
percentage  sand  figure,  secondly,  it  will  determine  the  corresponding 
numerical  values  for  these  three  soil  hydrological  parameters.  The 
values  are  then  automatically  stored  in  the  form  required  by  the 
infiltration  program  thus  reducing  the  amount  of  data  entry  required  of 
the  user.  The  program  to  generate  the  values  of  saturated  hydraulic 
conductivity  and  saturated  soil  moisture  content  has  been  developed  by 
the  SCS  at  Beltsville,  Maryland.  To  derive  the  saturated  hydraulic 
conductivity  for  example,  in  inches  per  hour,  the  following  expression 
is  used : 


2  2 

[-8. 9685-0. 0282( cl)+19.5235(POR)+0. 0001 (sd)  -0.0094(cl) 
e  2  2  2 

-8.3952(POR)  +0.0777(sd)(POR)-0.0029(sd)  (POR) 

2  2  2  2 
-0.0195(cl)  (POR)  -0.00002(sd)  (cl)-0.0273(cl)  (POR) 

2  2 
-0.0014(sd)  (POR)-0.000003(cl)  (sd)] 


Where : 


cl  -  percentage  clay 
sd  -  percentage  sand 
POR  -  porosity 
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The  Initial  moisture  content,  detention  capacity  and  iteration  period 
must  be  specified  by  the  user.  Again,  from  repeated  application  of  the 
model,  a  series  of  general  rules  will  be  derived  and  then  rather  than 
specifying  the  exact  numerical  figures  for  these  parameters,  the  user 
could,  by  supplying  a  more  general  level  of  information,  rely  on  the 
data  preparation  routines  in  the  model  to  derive  the  data  which,  on  the 
basis  of  past  experience,  are  considered  to  be  most  appropriate. 

Similarly,  the  precipitation  data  can  be  converted  to  the  format  which 
is  required  by  MILHY2,  from  the  form  in  which  they  are  available. 

Data  entry 

Under  the  proposed  scheme,  the  amount  of  data  entry  required  by  the  user 
is  reduced.  All  numerical  values  which  are  generated  by  the  data 
preparation  procedures  are  automatically  produced  in  the  form  required 
by  the  model. 

Data  checking 

It  is  necessary  to  check  the  data  before  model  execution  is  initiated. 

A  certain  degree  of  data  checking  can  also  be  incorporated  into  the 
program,  and  checks  on  units,  and  on  missing  or  incorrectly  typed  data 
will  certainly  be  very  effective. 

Figure  11  illustrates  the  nature  of  the  program  which  is  suggested  here. 
This  figure  illustrates  the  information  which  is  required  to  operate  the 
hydrograph  computation.  It  will  be  recalled  that  this  hydrological 
procedure  comprises  three  sections;  the  derivation  of  the  unit 
hydrograph,  the  derivation  of  incremental  runoff,  and  the  convolution  of 
these  two  series  to  produce  the  catchment  outflow  hydrograph.  Figure  11 
indicates  the  Information  which  must  be  supplied  by  the  user  and  the  two 
stages  of  data  preparation  and  checking  which  could  be  undertaken  by  the 
computer  program,  before  model  execution  begins.  Certainly  as  further 
enhancements  to  the  program  are  developed,  a  hierarchy  of  paths  through 
the  data  preparation,  entry  and*  checking  stages  could  be  provided 
depending  upon  the  nature  of  the  catchment  data  available,  and  the 
status  of  the  operator.  Further  refinement  could  involve  the 
incorporation  of  editing  facilities,  and  the  capability  to  view  and  to 
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In  this  series  of  applications  of  MILHY2  to  catchments  in  Vermont  and 
Iowa,  it  is  not  proposed  that  any  fine  tuning  of  the  model  parameters  be 
undertaken  to  assure. the  closest  fit  to  the  measured  hydrograph  which  is 
possible.  Rather,  the  catchment  data  which  have  been  derived  are  to  be 
used  in  one  application  to  each  storm.  Hence,  the  catchment  is  treated 
as  if  it  were  ungauged. 

To  assess  the  accuracy  of  the  model  predictions  for  this  wide  range  of 
experimental  frames,  the  same  two  stage  procedure  of  evaluation  will  be 
followed  (figure  2). 

In  total,  26  experimental  frames  (six  storms  applied  to  W-2,  North 
Danville,  Vermont  and  five  storms  to  each  of  the  four  catchments  in 
Treynor,  Iowa)  have  been  described  here.  Not  all  of  these  will  be 
reported  in  detail  during  the  following  discussion.  A  number  of 
selected  examples  will  serve  to  Illustrate  the  major  points  which  can  be 
made.  To  identify  each  experimental  frame,  the  catchment  name  and  the 
storm  number,  indicated  in  table  3,  will  be  provided. 

The  two  stage  procedure  which  compares  the  calculated  and  measured 
hydrographs  (figure  2)  will  be  followed  in  the  same  order  as  in  the 
comparison  of  the  predicted  hydrograhs  for  the  North  Creek  and  Sixmile 
Creek. 

Stage  I ;  Comparison  of  calculated  and  predicted  hydrograph 

A  comparison  of  calculated  and  measured  hydrographs  for  a  selection  of 
experimental  frames  is  provided  by  figures  12  to  16.  The  change  in 
scales  between  the  North  Danville  and  four  Treynor  catchments  should  be 
noted.  The  predictions  provided  by  MILHY2  for  W-2,  North  Danville  do 
not  approximate  the  measured  to  any  great  degree,  although  the  large 
vertical  scale  for  these  time  series  should  be  appreciated.  The  three 
storm  events  illustrated  in  figure  12  represent  the  range  of  inaccurate 
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Figure  12:  Comparison  of  calculated  and  measured  hydrographs  for  W-2 
North  Danville,  Vermont  (A)  Storm  3,  28  August  1970  (B) 
Storm  4,  16  July  1967  (C)  Storm  6,  2  June  1961 
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Figure  14:  Comparison  of  calculated  and  measured  hydrographs  for  W-2, 

Treynor,  Iowa  (A)  Storm  A,  20  June  1967  (B)  Storm  5,  7  June 
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and  inconsistent  results  which  are  obtained  for  this  catchment.  For 
storm  3,  (figure  12(A))  the  predicted  hydrograph  bears  no  similarity  in 
form  or  timing  to  the  measured.  Peak  discharge  is  also  highly 
overestimated.  The  measured  hydrograph  for  storm  4  (figure  12(B)) 
displays  a  double  peak.  The  calculated  hydrograph  also  has  a  double 
peak  but  neither  the  timing  nor  the  relative  magnitudes  of  the  two  peaks 
are  correct.  For  storm  6  (figure  12(C)),  the  model  predicts  a  much 
lower  runoff  than  was  experienced  in  the  catchment. 

MILHY2  provides  underpredictions  of  peak  discharge  for  all  10  storms 
applied  to  W-1  and  W-2,  Treynor,  and  figures  13  and  14  provide  four 
examples  of  this.  The  relationship  of  calculated  and  measured 
hydrographs  in  these  figures  is  very  similar  in  form  for  those  derived 
for  the  North  Creek  and  Sixmile  Creek  (DAJA37-81-C-0221 ) .  MILHY2  has  a 
tendency  to  overpredict  discharge  during  the  very  early  stages  of  the 
hydrograph  rise,  then  to  underpredict  discharge  during  the  peak  and 
finally  to  overpredict  discharge  during  the  latter  phases  of  recession. 
With  the  exception  of  storm  5  applied  to  W-1  however  (figure  13(B)),  the 
timing  of  the  predicted  hydrograph  quite  closely  approximates  the 
measured . 

Figure  15  provides  the  calculated  and  measued  hydrographs  for  storm 
numbers  3  and  4  applied  to  W-3  Treynor,  Iowa.  The  response  to  storm  3 
(figure  15(A))  is  typical  also  of  storms  1,  2  and  5  applied  to  this 
catchment.  The  measured  hydrograph  response  is  delayed  and  the  model 
does  not  predict  this.  The  overall  hydrograph  form  and  runoff  volume 
are  similar,  but  the  timing  is  poor.  The  prediction  for  storm  4  (figure 
15(B))  however  is  encouraging.  The  runoff  volume  and  timing  are  very 
well  predicted,  but  as  noted  above,  the  peaked  form  of  the  measured 
hydrograph  is  not  predicted  by  MILHY2.  Figure  16  illustrates  the 
overprediction  made  by  MILHY2  for  storm  4  on  W-4  Treynor,  Iowa  (figure 
16(A)).  The  predicted  response  to  storm  5  (figure  16(B))  again  has  a 
similar  relationship  to  the  measured  as  has  been  noted  for  the  North 
Creek  and  Sixmile  Creek. 

A  series  of  plots  of  calculated  against  measured  discharge  are  provided 
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by  figures  17  and  18.  The  dashed  line  indicates  the  position  of  perfect 
prediction  and  the  arrows  indicate  the  order  of  occurrence  of  errors 
from  t=0  and  at  successive  time  intervals  through  the  storm  event. 

Figure  17  illustrates  quite  clearly  the  range  of  overprediction  (storm 
3,  figure  17(A)))  to  underprediction  (storm  6,  figure  17(B))  derived  for 
this  catchment.  There  is  no  systematic  relationship  between  measured 
and  calculated  discharge  for  this  catchment.  The  patterns  of  hydrograph 
prediction  illustrated  in  figure  18(A)  for  storm  5,  W-1  and  in  figure 
18(B)  for  storm  5,  W-2,  Treynor,  Iowa  are  typical  of  the  response  to  the 
other  storms  applied  to  these  catchments,  and  are  also  similar  in  form 
to  those  produced  for  North  Creek  and  Sixmile  Creek  (figure  19).  A 

systematic  source  of  error  appears  to  occur  over  a  range  of  catchments 

which  causes  the  hydrograph  rising  limb,  peak  discharge,  and  beginning 
of  recession  to  be  underpredicted,  but  for  the  discharges  occurring 
during  the  latter  stages  of  recession  to  be  overpredicted. 

A  different  form  of  hydrograph  predictions  is  illustrated  for  storm  3 
applied  to  W-3  Treynor,  Iowa  in  figure  18(C).  Here,  the  pattern  is 
reversed,  overpredictions  of  the  rising  limb  and  underpredictions  of  the 
falling  limb  occur.  The  predicted  hydrograph  is  also  illustrated  to  be 
out  of  phase  with  the  calculated;  two  points  in  the  curve,  in  the  north 

and  east  corners,  are  observed  rather  than  the  more  usual  one,  in  the 

north  east  position.  Finally,  storm  5  applied  to  W-A  (figure  13(D)) 
displays  a  similar  pattern  to  the  Sixmile  Creek  and  North  Creek  where 
overprediction  of  the  rising  and  falling  limb  and  underprediction  of  the 
peak  discharge  have  produced  a  hydrograph  which  is  very  similar  in  terms 
of  runoff  volume,  but  not  as  peaked  as  the  measured. 

A  comparison  of  percentage  time  to  peak  discharge  error,  percentage  peak 
discharge  error,  and  percentage  mean  discharge  error  for  all  26 
experimental  frames  is  provided  in  figure  20.  Percentage  time  to  peak 
discharge  error  ranges  much  less  widely  than  the  other  two  indlcies. 

For  W-2,  North  Danville,  time  to  peak  discharge  is  predicted  exactly  for 
storm  A  and  underpredicted  for  the  other  five  storms  by  between  9%  and 
30%.  For  both  W-1  and  W-2,  Treynor,  the  exact  time  to  peak  discharge  is 
predicted  for  storms  2,  3,  and  A.  Storms  1  and  5  are  overpredicted  for 
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Relationship  between  discharge  predicted  by  HYM02  and 
measured  discharge  (A)  Storm  5,  7  June  1967,  W-1,  Treynor, 
(B)  Storm  5,  7  June  1967,  W-2,  Treynor  (C)  Storm  3,  14  June 
1967,  W-3,  Treynor  (D)  Storm  5,  7  June  1967,  W-4,  Treynor 
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Figure  19:  Relationship  between  discharge  predicted  by  HYM02  and  the 

measured  discharge  for  (A)  Storm  1,  9  October  1962,  North 
Creek  (B)  Storm  6,  4  May  1961,  Sixmile  Creek 


Figure  20:  Percentage  peak  discharge  error,  percentage  mean  discharge 

error,  and  percentage  time  to  peak  discharge  error  for  all 
26  experimental  frames 


both  catchments  by  between  9%  and  125%.  For  W-3  and  W-4,  percentage 
time  to  peak  discharge  error  ranges  from  -50%  to  +11%  and  -43%  to  +11% 
respectively.  Over  all  26  experimental  frames,  the  time  to  peak 
discharge  of  13  storms  are  predicted  to  within  plus  or  minus  10% 
(including  9  exactly)  and  only  in  4  cases  of  the  26,  is  the  prediction 
of  this  hydrograph  characteristic  in  error  by  greater  than  50%. 

Error  associated  with  peak  discharge  is  greater  than  that  for  time  to 
peak  discharge.  For  W-2,  North  Danville,  the  error  ranges  from  -82%  to 
+1882%  and  is  for  only  one  storm  within  20%  of  the  measured.  For  W-1 
and  W-2,  Treynor,  peak  discharge  is  underestimated  without  exception  by 
between  91%  and  67%.  For  W-3,  error  ranges  from  -43%  to  +689%. 

However,  the  greatest  range  of  error,  -33%  to  +3498%,  is  experienced  by 
W-4.  Over  all  26  experimental  frames,  there  are  no  events  where  peak 
discharge  is  predicted  to  within  10%.  In  fact,  in  19  of  the  26  cases, 
errors  of  greater  than  50%  occur. 

The  error  associated  with  the  prediction  of  mean  discharge  is  for  most 
storm  events  slightly  less  than  that  associated  with  peak  discharge. 
Very  wide  ranges  are  displayed  for  predictions  made  for  W-2,  North 
Danville,  and  W-3  and  W-4,  Treynor,  Over  all  26  experimental  frames, 
the  mean  discharge  of  three  storm  events  are  predicted  to  within  10% 
(including  two  exactly)  and  14  events  are  associated  with  error  of 
greater  than  50%. 

The  correlation  coefficients  and  error  standard  deviations  calculated 
for  these  26  experimental  frames  are  illustrated  in  figure  21.  The 
correlation  coefficients  are  very  low  and  indicate  very  little 
association  between  the  calculated  and  measured  hydrographs.  For  8  of 
the  26  cases,  a  correlation  coefficient  of  between  0.5  and  -0.2  exists, 
and  5  of  these  8  occur  for  W-2,  North  Danville.  Overall,  for  no  storm 
is  a  correlation  coefficient  of  greater  than  0.9  found.  The  error 
standard  deviation  values  Indicate  a  misleading  picture  of  better 
predictions  for  the  W-2  catchment.  North  Danville.  The  calculations  of 
this  statistic  are  affected  by  the  absolute  magnitude  of  the  discharges 
Involved,  and  which  for  this  catchment  are  indeed  very  small.  For  the 
Treynor  catchments  however  the  error  standard  deviations  are  still  low 


in  comparison  to  the  North  Creek  and  Sixmile  Creek,  a  maximum  of  2.7 
being  displayed. 

Stage  2:  Evaluation  of  errors 

Time  series  plots  of  model  forecast  error  (measured  discharge  minus 
calculated  for  each  time  interval)  for  a  selected  number  of  storms  are 
provided  in  figures  22  and  23,  for  each  catchment.  The  differences  in 
the  scales  of  the  vertical  axes  between  W-2,  North  Danville,  and  the 
Treynor  catchments  should  be  noted.  Much  less  error  is  associated  with 
the  prediction  of  the  small  discharges  measured  for  the  W-2,  North 
Danville  catchment. 

All  figures  confirm  the  tendency  (although  there  are  one  or  two 
exceptions)  towards  overprediction  (negative  error)  during  the  early 
stages  of  the  storm,  then  a  swing  upwards  to  underprediction  (positive 
error)  during  the  period  of  peak  discharge  and  a  tendency  back  to 
overprediction  during  the  latter  stages  of  recession.  A  similar  pattern 
in  errors  was  exhibited  for  the  North  Creek  (figure  24)  and  Sixmile 
Creek  (figure  25)  catchments. 

• 

A  plot  of  error  versus  the  measured  discharge  for  a  variety  of 
experimental  frames  is  provided  in  figure  26  for  W-2,  North  Danville  and 
in  figure  27  for  the  four  Treynor  catchments.  Figure  26  illustrates 
clearly  the  overprediction  for  storm  3  (figure  26(A))  and 
underprediction  for  storm  6  (figure  26(B)).  In  addition,  for  storm  6 
there  appears  to  be  an  almost  linear  relationship  between  error  and 
measured  discharge.  Indeed  these  two  series  have  a  correlation 
coefficient  of  0.99.  This  is  statistically  significant  at  the  95^ 
significance  level. 

In  figure  27,  all  four  plots  show  similar  systematic  forms  of  error  to 
the  North  Creek  and  Sixmile  Creek.  Storm  5  applied  to  W-1  (figure 
27(A))  and  W-2  (figure  27(B)). 

The  autocorrelation  functions  for  a  selection  of  representative  storms 
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for  each  catchment  are  indicated  in  figure  28.  All  of  these  functions 
indicate  a  much  lower  degree  of  autocorrelation  of  error  than  was  the 
case  for  the  North  Creek  and  Sixmlle  Creek.  Many  autocorrelation 
coefficients  approach  zero  by  lag  8.  However,  the  systematic  source  of 
error  in  model  prediction  is  still  significant. 

The  mean  and  standard  deviation  of  errors  is  provided  in  figure  29. 

Noticeably,  a  mean  very  close  to  zero  and  a  small  standard  deviation  are 

exhibited  by  North  Danville,  due  mostly  to  the  nature  of  the  small 

discharges  which  are  involved.  The  standard  deviation  of  error  is 

greatest  for  W-1  and  W-2,  where  one  standard  deviation  ranges  from  2.66 
3  -1 

to  0.8  m  s  .  For  W-3  and  W-4,  on  the  whole,  the  standard  deviations 

3  -1 

are  much  lower  (0.9  to  1.1  ms  ).  Over  all  26  experimental  frames,  17 

3  -1 

mean  errors  are  positive  and  range  from  0.1  to  1.44  m  s  indicating 

underprediction  by  the  model  (meausured  greater  than  calculated).  The 

3  -1 

negative  errors  range  from  -0.1  to  -1.08  m  s 

The  correlation  coefficients  in  table  4  indicate  that  for  none  of  the 
storms  documented  here  are  the  errors  normally  distributed. 

To  conclude  this  section  which  compared  the  predicted  and  measured 
hydrographs  for  a  variety  of  storms  and  for  5  catchments  in  Vermont  and 
Iowa,  the  following  two  points  can  be  made: 
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1  MILHY2  does  not  appear  to  provide  very  satisfactory  predictions  for 
W-2,  an  unnamed  tributary  of  the  Sleepers  River  catchment,  near  North 
Danville,  Vermont,  when  this  catchment  is  treated  as  an  ungauged 
catchment.  It  is  possible  that  improved  predictions  for  each  storm 
could  be  derived  if  a  degree  of  fine  tuning  of  the  model  parameters 
of  MILHY2  were  to  be  undertaken.  This  however,  is  not  the  point  of 
this  particular  exercise.  It  is  important  to  establish  the  degree  of 
accuracy  which  can  be  obtained  from  model  predictions  for  the 
ungauged  catchment.  Error  in  the  hydrograph  predictions  was  for  the 
North  Creek  and  Sixmlle  Creek,  attributed  to  model  and  data  error. 

The  likely  sources  of  model  error  in  the  context  of  the  application 
to  W-2,  North  Danville  will  now  be  examined. 
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Figure  29:  The  mean  (vertical  line)  and  one  standard  deviation 

(horizontal  bar)  of  discharge  error,  for  26  experimental 
frames 
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There  Is  a  large  probability  that  MILHY2  is  inappropriate  for 
application  to  this  particular  catchment.  Dunne  and  Black  (1970a, 
1970b)  document  observations  and  measurements  of  the  runoff  producing 
mechanisms  which  occur  in  a  small  area  of  the  Sleepers  River 
catchment  and  they  suggest  that  there  is  limited  evidence  to  suggest 
that  these  general  conclusions  may  be  extrapolated  for  most  of  the 
watershed.  The  major  runoff  producing  mechanism  is  overland  flow 
from  small  and  variable  contributing  areas  located  adjacent  to  the 
stream,  in  poorly  drained  positions  where  the  water  table  is  near  to 
the  surface.  Runoff  from  these  areas  reaches  the  channel  very 
quickly.  MILHY2  is  not  designed  to  model  these  particular 
hydrological  processes  in  terms  of  the  methods  used  to  generate 
runoff  and  the  use  of  unit  hydrograph  procedures  to  route  this  runoff 
through  the  catchment  area.  Hortonian  overland  flow  occuring  over 
large  areas  has  not  been  observed  on  this  catchment  and  indeed,  the 
Infiltration  capacity  of  the  soils  exceeds  most  measured  rainfall 
intensities . 

There  is  not  such  a  high  probability  that  data  errors  will  be  large 
for  this  catchment.  As  an  ARS  experimental  watershed,  it  is  likely 
that  precipitation  and  measured  hydrograph  information  will  be  as 
reliable  as  possible.  It  is  possible  however,  that  the  soils  data 
which  are  derived  from  the  Brakensiek  and  Rawls  charts  are  not 
accurate  for  simulation  in  this  small  catchment. 

2  For  the  four  catchments  located  near  to  Treynor,  Iowa,  again  when 
they  are  treated  as  ungauged  catchments,  a  wide  range  of  predictions 
is  derived.  Overall,  very  similar  patterns  (but  not  magnitude)  of 
discharge  prediction  error  are  obtained  as  were  derived  from 
application  to  the  North  Creek  and  Sixmile  Creek.  The  timing  of  the 
predicted  hydrographs  is  good,  but  peak  discharge  is  commonly 
underpredicted  and  a  systematic  source  of  error  is  identified,  where 
mean  errors  differ  from  zero,  are  not  normally  distributed,  and 
exhibit  autocorrelation. 

Again,  improvements  to  the  unit  hydrograph,  the  most  likely  source  of 


such  systematic  error,  can  be  suggested.  Certainly,  the 
dimensionless  unit  hydrograph  method  which  is  used  by  MILHY2  has  not 
been  calibrated  for  catchments  containing  contour  corn,  located  in 
Iowa,  whereas  it  has  been  for  Texas  and  Arkansas.  This  feature  may 
also  be  connected  with  the  scale  of  the  catchments.  It  is  possible 
that  better  predictions  will  be  derived  for  larger  catchments  than 
the  small  ones. 
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Figure  3 1 : 


A  comparison  of  Che  change  in  moisture  context  at  6  minute 
intervals  which  are  predicted  by  the  infiltration  model  for 
the  Ida  silt  loam,  and  associated  with  the  application  of  a 
storm  of  22  June  1964  (total  precipitation  27.94  mm)  for 
(A)  a  30  second  iteration  period  (B)  a  10  second  iteration 
period  (C)  a  10  second  iteration  period  and  halved  cell 
dimensions 


Few  cases  of  physically  unrealistic  infiltration  behaviour  were 
experienced  in  any  application  of  MILHY2  which  has  been  considered  in 
this  report.  Unrealistic  behaviour  can  be  demonstrated  to  occur  in 
association  with  a  combination  of  very  small  cell  size  in  the  soil 
column,  small  time  Increments,  and  high  precipitation  intensity. 

Figure  30  illustrates  the  precipitation  and  resulting  infiltration  and 
runoff  behaviour  for  all  five  soil  types  in  the  W-2,  North  Daville, 
Vermont  for  storm  number  4.  Infiltration  is  represented  by  the  changing 
moisture  content  of  the  five  soil  columns  at  three  depths,  0.05  metres, 
0.15  metres  and  0.3  metres  every  30  minutes  from  04:30  hours  (the  start 
of  the  storm),  for  9  hours  (storm  duration).  For  each  soil  type,  the 
most  rapid  and  greatest  increase  in  soil  moisture  content  is  experienced 
at  the  shallowest  depth  indicated.  The  increase  in  soil  moisture 
content  further  down  the  soil  column  is  not  as  great,  and  occurs  more 
gradually.  Runoff  occurs  in  association  with  saturated  surface 
conditions  and  higher  rainfall  intensities.  Where  a  greater  amount  of 
precipitation  is  required  to  saturate  the  soil  (Colrain  compared  to 
Peacham,  for  example),  less  runoff  results. 

Figure  31  illustrates  the  effect  which  the  choice  of  the  cell  size  and 
iteration  period  has  upon  infiltration  behaviour,  again  as  represented 
by  changes  in  soil  moisture  content.  These  results  were  derived  from 
application  of  a  storm  of  22  June,  1964  (which  has  not  previously  been 
used  in  this  thesis)  which  has  a  total  of  27.94  mm  precipitation  to  the 
soil  column  Ida  (a  silt  loam)  which  occurs  in  the  watersheds  near 
Treynor,  Iowa.  This  soil,  in  the  absence  of  more  detailed  data,  is 
assumed  not  to  be  layered  and  is  represented  by  a  soil  column  comprising 
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6  cells.  The  hydrological  characteristics  have  been  derived  from  the 
centroid  position  on  the  Brakensiek  and  Rawls  charts.  Figures  31(A), 
31(B),  and  31(C)  all  illustrate  the  initial  moisture  content  and  the 
moisture  content  at  successive  6  minute  intervals  for  each  cell.  Figure 
31(A)  Illustrates  the  response  when  a  30  second  iteration  period  is 
assumed;  figure  31(B)  if  a  10  second  period  is  assumed;  and  figure  31(C) 
where  both  a  10  second  iteration  period  and  twice  as  many  cells,  with 
halved  cell  dimension  are  used.  There  is  very  little  difference  between 
the  soil  moisture  content  profiles  which  develop  during  the  storm  when 
the  6  cells  are  utilized,  and  iterations  of  30  or  10  seconds  are  used. 
Halving  the  cell  size,  however,  has  no  effect  during  the  first  4  time 
intervals,  but  during  the  next  3  time  intervals,  a  form  of  physical 
instability  occurs  and  moisture  content  oscillates  through  a 
range  of  0.2  m  m  .  This  instability  corresponds  to  periods  where  large 
amounts  of  precipitation  occur.  When  the  precipitation  amount  drops 
again,  for  intervals  8  to  10,  the  profile  resumes  a  physically  realistic 
form  and  one  which  is  similar  to  those  attained  in  figures  31(A)  and 
31(B).  It  is  interesting  to  note  that  associated  with  these  conditions 
is  a  value  of  (BAL)  (equation  2),  a  measure  of  the  mean  numerical  error, 
of  0.015  for  condition  'C'  compared  with  a  value  of  0.010  for  condition 
'A'.  No  benefit  is  seen  to  be  derived  from  the  adoption  of  smaller  cell 
sizes  and  shorter  time  increments. 


BAL  =  0  -  0  -  ci  +  ce  +  cd 

end  init 


(2) 


Where : 

3  -3 

BAL  -  numerical  error  (ram  ) 

3  -3 

0  -  total  water  content  of  soil  profile  (mm  )  at  end  of 

end  3  -3 

0  -  initial  total  water  content  of  entire  profile  (mm  ) 

init 

ci  -  cumulative  infiltration  (ms  ) 
ce  -  cumulative  evaporation  (ms  ) 
cd  -  cumulative  drainage  (ms  ) 


Slightly  higher  errors  are  exhibited  for  more  complex  soil  and 
precipitation  conditions.  Table  5  provides  the  details  of  the  value  of 
(BAL)  (a  measure  of  the  magnitude  of  numerical  errors  incurred  by  the 
solution  of  the  Richards  equation  using  an  explicit  finite  difference 
method)  for  each  soil  type  on  all  seven  catchments  located  in  Texas, 
Arkansas,  Vermont  and  Iowa  for  all  storms  which  have  now  been 
documented.  For  many  cases,  the  value  of  (BAL)  can  be  related  to  soil 
depth,  soil  type,  and  precipitation  intensity.  For  example,  the  results 
presented  in  table  5  for  North  Creek,  Texas  illustrate  that  greater 
errors  occur  for  the  soil  column  representative  of  the  Gowen-Pulexas 
soil  groups.  This  soil  column  is  deeper  than  those  representing  the 
Bonti-Cona-Truce  and  Thurber-Hasse  soil  groups,  and  consequently  has  a 
greater  number  of  cells  for  which  a  solution  must  be  provided.  The 
Gowen-Pulexas  also  has  a  higher  conductivity  than  the  other  two  soils, 
which  both  have  clay  in  layers  2  and  3  (tables  6,  7,  and  8).  The  lowest 
error  for  the  Gowen-Pulexas  soil  occurs  for  storm  3.  This  storm  has  the 
shortest  duration  (1.3  hours)  and  the  most  precipitation  (107  mm).  In 
contrast,  the  greatest  error  for  this  soil  type  occurs  for  storm  1  which 
is  8.25  hours  long  and  throughout  is  very  erratic;  periods  of  high 
precipitation  intensity  alternate  with  periods  of  very  little  rain. 

Such  rapid  fluctations  in  rainfall  intensity  in  successive  time 
intervals  appear  to  be  associated  with  greater  errors  in  the  solution  of 
the  Richards  equation. 

Very  similar  relationships  between  soil  characteristics  and  the  value  of 
(BAL)  are  exhibited  by  the  information  provided  for  the  storms  applied 
to  the  Sixmile  Creek.  Larger  errors  are  associated  with  the  deeper 
soil,  Leadvale.  However,  for  this  suite  of  storms,  there  is  no  clear 
relationship  between  (BAL)  and  storm  characteristics. 

For  W-2,  North  Danville,  the  magnitude  of  error  is  very  much  less  than 
has  been  noted  for  the  previous  two  catchments.  This  may  be  related  to 
the  shallow  soil  columns  which  were  used  to  represent  the  soils  of  this 
catchment.  The  greater  amount  of  numerical  error  is  not  consistently 
associated  with  the  same  soil  column.  The  Cabot  soil  type  exhibits  the 
greatest  error  for  storms  1,  4,  and  5,  and  the  Woodstock  soil  type  for 


•/‘.V.V 


P  m 

.  ^  »  •  .  '  . 


V  -•  ■ 


North  Creek,  Texas 

Gowen-Pulexas  Bonti-Cona-Truce  Thurber-Hasse 


Slxmile  Creek,  Arkansas 

Leadvale 

Enders 

Mountainbur 

1  -0.6 

-0.2 

-0.2 

2  -0.7 

-0.2 

-0.2 

3  -3.6 

-0.4 

-1.0 

4  -4.3 

-0.2 

-0.1 

5  -1.1 

-0.6 

-0.5 

6  -0.6 

-0.1 

-0.8 

W-2;  North  Danville, 

Vermont 

Colrain  Peacham  Calais  Cabot  Woodstock 

1  0.0  -0. 

1  0 

.0 

1.1 

-0.5 

2  0.0  0. 

0  0 

.0 

0.0 

-0.4 

3  0.0  0. 

0  0 

.0 

0.3 

-0.5 

4  -0.2  0. 

0  0 

.0 

1.3 

-1.2 

5  0.0  -0. 

3  -0 

.1 

1.6 

0.0 

6  0.0  0. 

0  0 

.0 

0.0 

-0.2 

W-1,  Treynor,  Iowa 


Marshall 
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Table  5  ...  continued  from  previous  page 


*  -2^-3 

BAL  (xlO  mm  ) 

'-1*  Storm  _ 

number 

Soli  types 

i  - ^ - 


W-2. 

Mona 

Treynor,  Iowa 

Marshall 

Napier 

Ida 

1 

-0.3 

0.0 

0.0 

-2.3 

2 

-0.1 

0.0 

0.0 

-0.9 

3 

0.0 

0.0 

0.0 

-0.8 

4 

-3.3 

0.0 

0.0 

-11.9 

5 

-0.3 

0.0 

0.0 

-3.0 

W-3,  Treynor,  Iowa 


1 

-0.3 

0.0 

0.0 

-2.8 

2 

0.0 

0.0 

0.0 

-1.4 

3 

0.0 

0.0 

0.0 

-0.9 

•*  _ 

4 

-2.0 

0.0 

0.0 

-9.2 

•*. 

5 

-0.1 

0.0 

0.0 

-0.8 

W-4, 

Treynor,  Iowa 

k'- 

1 

-2.7 

0.0 

0.0 

-2.8 

2 

0.0 

0.0 

0.0 

-0.9 

f  ' 

3 

0.0 

0.0 

0.0 

-0.9 

m 

4 

-2.0 

0.0 

0.0 

-9.2 

5 

-0.1 

0.0 

0.0 

-0.8 

f  ' 

m 


the  remaining  three  storms.  These  two  soils  do  not  have  any  particular 
characteristics  in  common,  and  the  deepest  soil  for  this  catchment  with 
the  greatest  number  of  cells  is  Colrain. 

For  all  4  catchments  near  Treynor,  the  soil  column  representing  the  Ida 
soil  type  exhibits  the  greatest  error.  This  soil  column  is  the 
shallowest,  but  the  cell  dimensions  are  the  smallest.  For  all  four 
catchments,  the  greatest  error  is  experienced  for  storm  4.  This  storm 
has  the  highest  precipitation  total,  but  also,  as  noted  for  Jexas,  the 
most  rapidly  alternating  successions  of  high  and  low  intensity  rainfall. 
The  lowest  error  for  W-1  and  W-2  is  associated  with  storm  3  which  has 
the  lowest  total  precipitation.  The  lowest  error  for  W-3  and  W-4  is 
associated  with  storm  5  which  has  the  second  lowest  pecipitation  total, 
but  the  shortest  duration. 

The  relationship  of  error  to  precipitation  is  demonstrated  in  figure  32. 
The  information  for  this  figure  is  taken  from  storm  4  applied  to  W-1, 
Treynor.  Cumulative  precipitation  is  compared  to  cumulative  (BAL)  for 
the  two  soil  columns  which,  as  indicated  in  table  5,  exhibit  errors  in 
solution.  A  steeper  gradient  on  the  cumulative  precipitation  curve 
appears  to  be  related  to  a  steeper  rise  in  the  value  of  cumulative  BAL 
for  each  soil  type.  Indeed,  the  correlation  coefficient  between 
cumulative  precipitation  and  the  cumulative  (BAL)  for  Monona  soil  type 
is  0.964  and  for  the  Ida  soil,  is  0.997.  Both  of  these  correlation 
coefficients  are  slgnficant  at  the  95%  confidence  level. 

Over  all  experimental  frames,  it  is  not  considered  that  numerical  errors 
are  large  enough  to  justify  an  examination  of  alternative  numerical 
techniques . 


Summary  Of  Applications 


4.1  Introduction 


To  summarize  the  results  of  the  application  of  MILHY2  to  38  storms,  and 
for  a  range  of  seven  catchments  in  Texas,  Arkansas  (2),  Vermont,  Iowa, 
figures  33,  34,  and  35  have  been  produced.  Figure  33  attempts  to  assess 
the  accuracy  of  MILHY2  for  the  prediction  of  peak  discharge;  figure  34, 
the  accuracy  of  the  time  to  peak  discharge  predictions  and  figure  35, 
the  closeness  of  the  overall  hydrograph  form.  From  these  figures,  the 
following  conunents  may  be  derived: 


4.2  Prediction  of  peak  discharge 


Figure  33(A)  provides  a  plot  of  calculated  versus  measured  peak 
discharges  for  all  38  experimental  frames.  A  correlation  coefficient  of 
0.911  between  these  two  series  has  been  calculated.  This  is  not 
statistically  significant,  and  the  trend  towards  underprediction  of  peak 
discharge,  which  has  been  noted  previously,  is  seen  clearly.  This  type 
of  plot,  although  often  produced  in  modelling  studies,  is  slightly 
misleading  in  that  the  very  small  deviations  from  the  dashed  line 
(indicating  perfect  prediction)  in  the  lower  peak  discharge  range  can 
be,  in  relative  terms,  a  good  deal  more  significant  than  the  apparently 
larger  deviations  which  occur  at  higher  discharges.  This  point  is 
lll.ustrated  by  figure  33(B),  where  percentage  peak  discharge  error 
plotted  against  measured  discharge  is  given  by: 
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Figure  33;  a  suonary  of  Che  accuracy  of  HYM02  for  .the  prediction  of 
peak  di^harge  over  all  38  experimental  frames  (A)  the 
relationship  of  calculated  and  measured  peak  discharge  (6) 
Che  relationship  of  percentage  peak  error  and  measured  peak 
discharge  (C)  the  relationship  of  percentage  peak  discharge 
error  and  total  precipitation 
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Figure  34: 


Total  pr$!cipitation  fm  m) 


A  summary  of  the  accuracy  of  HYM02  for  Che  prediction  of 
the  time  to  peak  discharge  over  a^l  38  experimental  frames 
(A)  the  relationship  of  calculated  and  measured  time  to 
peak  discharge  (B)  the  relationship  of  percentage  time  to 
peak  discharge  error  and  measured  peak  discharge  (C)  the 
relationship  of  percentage  time  to  peak  discharge  error  and 
total  precipitation 
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PDE  =  q  m  -  q  c  x  100%  (3) 

P  P 


q  m 
P 


Where : 


q  m 
P 

q  c 

P 

Much  greater  error  Is  seen  to  be  associated  with  the  prediction  of  lower 
peak  discharge  than  with  higher.  Indeed,  this  figure  suggests  that  the 
closest  estimate  of  peak  discharge,  provided  by  MILHY2,  will  be  derived 
for  peak  discharges  between  the  range  20  to  65  m  s  .  There  is  a 
greater  tendency  towards  overestimation  within  the  lower  discharges,  and 
underestimation  at  higher. 

Figure  33(C)  provides  a  plot  of  percentage  peak  discharge  error  versus 
total  precipitation.  From  this  range  of  experimental  frames,  there  does 
not  appear  to  be  a  clear  relationship  between  these  two  series.  However, 
it  could  be  suggested  that  in  general,  greater  accuracy  is  provided  by 
MILHy2  for  the  prediction  of  the  peak  discharge  for  larger  storms. 

4.3  Predictions  of  time  to  peak  discharge 


-  measured  peak  discharge  (ft  s  ) 

3  -1 

-  calculated  peak  discharge  (ft  s  ) 


MILHY2  predicts  the  time  to  peak  discharge  much  more  accurately  than  any 

other  hydrograph  characteristic.  The  correlation  between  calculated  and 

measured  time  to  peak  discharge,  indicated  in  figure  34(A),  is  0.974. 

This  is  higher  than  that  calculated  for  the  association  between 

calculated  and  measured  peak  discharge.  Figure  34(B)  indicates  that 

over  the  total  range  of  measured  peak  discharges  which  are  considered  in 

this  study,  a  much  lower  percentage  error  for  time  to  peak  discharge  is 

derived,  than  for  peak  discharge.  There  are  just  one  or  two  outliers, 

3  -1 

for  example  at  12  m  s  .  This  can  be  identified  as  the  error  associated 
with  the  prediction  of  time  to  peak  discharge  for  storm  4,  W-1 ,  Treynor. 
As  the  other  errors  for  this  hydrograph  characteristic  are  much  lower. 


this  outlier  might  possibly  be  associated  with  error  in  the 
precipitation  or  measured  hydrograph  data  which  were  utilized  for  this 
particular  storm  event.  Figure  34(C)  also  indicates  very  little  clear 
relationship  of  percentage  time  to  peak  discharge  error  to  precipitation 
totals. 

4.4  Predictions  of  the  overall  form  of  the  discharge  hydrograph 


The  closeness  of  form  of  the  calculated  to  measured  hydrograph  is,  for 

the  purposes  of  this  comparison,  indicated  by  the  value  of  the 

correlation  coefficient.  Figure  35(A)  provides  the  distribution  of  the 

correlation  coefficient  according  to  measured  peak  discharge.  On  the 

whole,  a  closer  association  is  derived  for  hydrograph  events  where  peak 

3-1 

discharge  ranges  between  20  and  60  m  s  .  Below  and  above  these  values, 
the  correlation  coefficient  between  the  calculated  to  measured  increases 
in  range.  Figure  35(B)  indicates  no  clear  relationship  between  the 
correlation  coefficient  and  total  storm  precipitation,  although  very 
generally,  the  closeness  of  fit  does  have  a  tendency  to  improve  as  the 
total  precipitation  increases. 

MILHY2  does  also  appear  to  provide  more  accurate  predictions  for  some 
catchments  than  others.  To  assess  the  overall  goodness  of  fit  of  the 
calculated  hydrographs  for  the  range  of  storms  applied  to  each 
catchment,  a  multiple  index  (I  )  was  derived  from  the  percentage  peak 

X 

discharge  error  (PDE),  percentage  time  to  peak  error  (TPE) ,  and  the 
correlation  coefficient  (r)  according  to  the  following  expression: 


PDE  I  +  I  TPE  I  +  lOO(l-r) 


This  index  was  evaluated  for  each  experimental  frame,  and  the  mean  value 
was  derived  for  each  catchment.  The  results  of  this  are  presented  in 
table  9.  For  the  range  of  storms  which  have  been  considered  in  this 
analysis,  the  best  predictions  are  derived  for  the  Sixmlle  Creek, 
Arkansas,  and  then  for  the  North  Creek,  Texas.  The  model  does  not 
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Catchment 

Value  of 

M 

< 

O 

each 

storm 

Mean 

A 

value 

1 

2 

3 

4 

5 

6 

of  I 

X 

North  Creek, 

Texas 

62 

45 

150 

104 

9 

42 

69 

Sixmlle  Creek, 
Arkansas 

62 

18 

7 

24 

27 

23 

27 

W-2,  North 

Danville,  Vermont 

69 

402 

11961 

71 

139 

211 

2142 

W-1,  Treynor, 

Iowa 

196 

149 

139 

229 

117 

166 

W-2,  Treynor, 

Iowa 

188 

104 

117 

125 

116 

130 

'W-3,  Treynor, 

Iowa 

758 

185 

159 

69 

47 

244 

W-A,  Treynor, 

Iowa 

3557 

28 

80 

322 

55 

808 

appear  to  provide  suitable  predictions  for  the  unnamed  tributary,  W-2, 
of  the  Sleepers  River  catchment.  In  comparison  to  this  catchment,  it 
was  more  successful  for  the  four  catchments  near  Treynor,  Iowa.  In  this 
context,  it  should  be  recalled  that  the  unit  hydrograph  procedure  has 
been  calibrated  for  34  catchments  located  in  Texas,  Oklahoma,  Arkansas, 


Louisiana,  Mississippi,  and  Tennessee 


Application  of  MILHY2  has  provided  a  range  of  results  from 
catchments  in  Texas  and  Arkansas  (see  (2)),  and  Vermont  and  Iowa 
(this  report).  The  following  points  are  worthy  of  note: 

i)  the  correlation  between  predicted  and  measured  peak 
discharge  using  MILHY2  Is  high  (r  =  0.91) 

11)  the  time  to  peak  dsicharge  estimation  is  particularly  good 
using  MILHY2  (correlation  between  predicted  and  measured 
-  0.97) 

ill)  the  prediction  for  w-2  (Sleepes  River  Catchment)  is  poor 
(see  figure  12) 

Iv)  comparison  of  MILHY2  and  MILHY  (HYMO)  f 'r  32  experimental 
frames  shows  strong  evidence  of  the  overall  improvements 
achieved  by  MILHY2  (figures  36  and  37),  especially  in  time 
to  peak  discharge 

It  is  recommended  that  further  field  trials  of  MILHY2  are 

undertaken  (this  work  is  currently  taking  place  under 
DAJA45-85-C-0022)  and  that  the  computing  needs  of  MILHY2  are 
explored  with  respect  to  run-time  performance. 
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36:  Comparison  of  the  percentage  time  to  peak  discharge  error 

ofMILHY2andHYMO,  for  32  experimental  frames 
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Figure  37:  Comparison  of  the  percentage  peak  discharge  error  of  IVIILHY2 

and  HYMO,  for  32  experimental  frames 


non  no  noon 


Program:  MILHY2 

HYMO  including  a  physically  based  infiltration  algorithmn 
which  replaces  the  Soil  Conservation  Service  curve  number 
model 

Coded  by:  S  Howes 

University  of  Bristol 

Notes:  Much  of  the  code  remains  unaltered  but  a  number  of 

subroutines  and  functions  have  been  added. 

All  additional  code  is  written  in  FORTRAN?? 


C 

C 

C 


Modifications  occur  in  following  subroutines: 
CMHYD 
ERROR 


C 

C 

C 

C 

C 

C 

C 


Additional  subroutines: 
SOILM 
HYDCON 
TWO 
GRAD 
SMCURV 
BLOCK  DATA 


C  Additional  functions: 

C  RMAX 

C  RMIN 


C 


OPEN  (1 ,STATUS="OLD",FORM="FORMATTED",FILE="datal",MODE="IN") 
OPEN(25,FORM="FORWATTED",FILE="data2",MODE="IN",STATUS="OLD") 
0PEN(6,F0RM=''F0RMATTED",STATUS="NEW",M0DE=”0UT”,FILE==’’results") 

COMMON /BLOCKl/  OCFS(300,6) ,DATA(310) ,CFS(300) ,CTBLE(50 , 1 1) , 
&RAIN(300) ,R0IN(6) , 

&A(20,6),Q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 
&ZALFA(20),IEND(6:i,DA(6),DIST(6),SEGN(6),DT(6),PEAK(6),ISG(6), 
&NPU , NHD , NER ,MAXNC , NCOMM , ICC , NCODE ,TIME , KCODE , ICODE 


C  Definition  of  variables  in  common 


C  OCFS 
C  DATA 
C  CFS 
C  CTBLE 
C  RAIN 
C  ROIN 
C  A 
C  Q 

C  DEEP 
C  ITBLE 
C  DP 


Hydrograph  discharge 
Data  associated  with  each  command 
Unit  hydrograph  discharge 
Command  table 

Cumulative  precipitation  values 
Runoff  volume  of  discharge  hydrograph 
End  area 

Flow  rate  for  rating  curve 

Elevation  of  water  surface  (for  rating  curve) 

Integer  table 

Flow  depth  for  previously  computed  travel  time  flow  relationship 


C  SCFS  Discharge  for  previously  computed  travel  time  flow  relationship 
C  C  Travel  time  coefficient  for  previously  computed  travel  time 

C  flow  relationship 

C  ZALFA  Alphnuraeric  code  table 
C  lEND  Number  of  points  in  the  hydrograph 

C  DA  Drainage  area 

C  DIST  Segment  boundary  point  for  each  segment  of  a  cross  section 
C  SEGN  Mannings  'n'  for  each  segment  of  a  cross  section 
C  DT  Time  increment  for  rainfall  or  discharge 

C  PEAK  Peak  discharge  for  hydrograph 

C  ISG 

C  NPU  Punch  code 

C  NHD  Hydrograph  identification  number 

C  NER  Error  number 

C  MAXNO  Maximum  number  of  data  entires  to  be  expected  for  any  command 
C  NCOMM  ,  Number  of  commands 
C  ICC  Continuation  line 

C  NCODE  'Number  of  command 

C  TIME  Start  time  of  simulation 

C  KCODE  Measurement  unit  of  input 

C  0  -  imperial 

C  not  0  -  metric 

C  ICODE  Measurement  unit  of  output 
C  0  -  imperial 

C  not  0  -  metric 

NCODE=0 

NPU=0 

ICC=0 

1  NER-0 
CALL  HONDO 

IF  (NER)  2,2,19 

2  GO  TO  (3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19),  NCODE 

3  T1ME=DATA( 1 ) 

NPU=DATA(2) 

KCODE=DATA(3) 

IC0DE=DATA(4) 

GO  TO  1 

4  CALL  STHYD 

GO  TO  1 

5  CALL  RECHD 

GO  TO  1 

6  CALL  CMPHYD 

GO  TO  1 

7  CALL  PRTHYD 

GO  TO  1 

8  CALL  PUHYD 

GO  TO  1 

9  CALL  HPLOT 

GO  TO  1 

10  CALL  ADHYD 

GO  TO  1 

1 1  CALL  SRC 
GO  TO  1 

12  CALL  CMPRC 

GO  TO  1 
CALL  STT 


13 
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14 

GO  TO  1 
CALL  CMPTT 

15 

GO  TO  1 
CALL  ROUTE 

16 

GO  TO  1 
CALL  RES VO 

17 

GO  TO  1 
CALL  ERROR 

18 

GO  TO  1 
CALL  SEDT 

19 

GO  TO  1 
STOP 

END 

SUBROUTINE  HONDO 

C  This  subroutine  reads  in  the  data  from  'datal',  searches  an  alphanumeric 
C  code  table  to  determine  the  NCODE  of  the  required  operation,  and  collects 
C  variables  from  the  freefloating  data  field. 

C  The  command  table  (CTBLE),  Integer  table  (ITBLE),  number  of  commands 
C  (NCOMM)  and  alphanumeric  array  (ZALFA)  are  initialized  in  BLOCK  DATA 
C  located  at  the  end  of  this  listing. 


COMMON/BLOCKl /  OCFS ( 300 , 6 ) , DATA(3 10 ) ,CFS (300 ) , CTBLE ( 50 , 1 1 ) , 
&RAIN(300),ROIN(6), 

&A(20,6),Q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 

&ZALFA(20),IEND(6),DA(6),DIST(6),SEGN(6),DT(6),PEAK(6),ISG(6), 

&NPU , NHD , NER , MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 

DIMENSION  CHAR(60),  ALPHA( 1 1 ) ,AUXA( 10) ,AUXB( 10) 

IF  (ICC)  1,1,3 
C  READ  IN  DATA  CARD 

1  READ  (1,42)  (ALPHA(I),I=l,ll),(CHAR(I),I=l,60) 

IF  FIRST  CHARACTER  IS  BLANK  THE  CARD  IS  A  CONTINUATION  OF 
PREVIOUS  CARD. 

IF  (ALPHA(l)-ZALFA(ll))  2,9,2 
IF  (ICC)  3,3,40 

ASTERISK  IN  COL.  80  MEANS  SKIP  TO  NEW  PAGE  BEFORE  PRINTING  CARD 
IF  (CHAR(60)-ZALFA(11))  4,5,4 
WRITE  (6,43) 

WRITE  (6,44)  (ALPHA(I),I=1,11),(CHAR(I),I=1,60) 

IF  FIRST  CHARACTER  IS  A  *  THE  PREVIOUS  CARD  WAS  A  COMMENT  CARD 
IF  (ALPHA(1)-ZALFA(12))  10,6,10 

IF  PUNCH  CODE  POSITIVE,  COMMENT  CARDS  ARE  PUNCHED. 

IF  (NPU)  8,8,7 

WRITE  (7,45)  (ALPHA(I),I=1,11),(CHAR(I),I=1,60) 

ICC-O 
GO  TO  1 

9  WRITE  (6,44)  ( ALPHA( I ) , 1=1 , 11 ) , (CHAR( I ) , 1=1 , 60) 

GO  TO  24 

C  SEARCH  FIRST  TWO  ALPHAMERIC  CHARACTERS  TO  SEE  IF  THEY  ARE  NUMBERS 


DO  12  I=«l,10 

IF  (ALPHA(1)-2ALFA(I))  11,15,11 
IF  (ALPHA(2)-2ALFA(I))  12,15,12 
CONTINUE 

STATEMENT  NUMBER  7  IS  BRANCHED  TO  IF  NUMBERS  ARE  PRESENT 
IF  NOT  NUMBER  SEARCH  COMMAND  TABLE  FOR  MATCH 
CALL  FIRST  10  VALUES  FROM  PERMANENT  DATA  STORAGE 
DO  14  I-l,NCOMM 
DO  13  J=l,ll 

IF  (CTBLE(I,J)-ALPHA(J))  14,13,14 

SN  10=PART  MATCH 

CONTINUE 

IF  THIS  LOOP  IS  COMPLETED  WE  HAVE  COMPLETE  MATCH-  CALL  NCODE 

AND  MAX  NUMBER  AND  EXIT  LOOP 

NCODE=ITBLE(I,l) 

MAXNO=ITBLE(I,2) 

GO  TO  21 
CONTINUE 

IF  MAJOR  LOOPS  FINISHED  WITHOUT  A  MATCH  WRITE  ERROR  MESSAGE 

AND  SET  NER  =  1 

NER=1 

WRITE  (6,46) 

RETURN 

CONVERT  DIGIT  INPUT  CODE  FROM  ALPHAMERIC  TO  INTEGER  FORM 

NCODE=GIT (ALPHA, 1,2,1. )+0. 5 

FIND  MAX  NUMBER  OF  DATA  ITEMS  FOR  THIS  NCODE 

DO  17  I=1,NC0MM 

IF  (ITBLEd,  1)-NC0DE)  17,16,17 

MAXNO-ITBL£(r,2) 

GO  TO  21 
CONTINUE 

SEARCH  DATA  ROUTINE 

SEE  IF  ANY  DATA  FOR  THIS  CARD 

DO  19  I=1,NC0MM 

IF  (ITBLEd,  l)-NCODE)  19,18,19 

MAXN0=ITBLE(I,2) 

GO  TO  20 
CONTINUE 
CONTINUE 

IF  (MAXNO)  23,22,23 
RETURN 

ZERO  ARRAYS  AND  COUNTERS 
DO  47  1=1,310 
DATA  (I)=0. 

NDATA=1 
NCHAR=0 
DO  26  1-1,10 
AUXA(I)=0. 

AUXB(I)=0. 

IT1  =  1 
IT2-1 
SIGN-1. 

LDGIT-0 

KDGIT=0 

CARRY  OUT  DIGIT  BY  DIGIT  SEARCH  AND  ACCUMULATION 
NCHAR-NCHAR+1 

HAVE  WE  CONSIDERED  ALL  CHARACTERS  -  RETURN  IF  SO 


IF  (NCHAR-60)  28,32,1 
DO  29  1=1,15 

IF  (CHAR(NCHAR)-ZALFA(I))  29,30,29 

CONTINUE 

GO  TO  32 

GO  TO  (33,33,33,33,33,33,33,33,33,33,32,27,36,32,31,27),  I 
SN  39  HANDLES  SIGN  CONTROL  ON  1130  VERSION 
SIGN=-1.0 
GO  TO  27 

CHARACTER  IS  BLANK  OR  COMMA  -  DOES  IT  FOLLOW  A  DIGIT 
GO  TO  (27, A8),  ITl 

CHARACTER  IS  A  DIGIT  -  HAS  A  DECIMAL  BEEN  ENCOUNTERED 

GO  TO  (34,35),  IT2 

LDGIT=LDGIT+1 

IT1=2 

AUXA(LDGIT)=CHAR(NCHAR) 

GO  TO  27 

KDGIT=KDGIT+1 

AUXB(KDGIT)=CHAR(NCHAR) 

GO  TO  27 

CHARACTER  IS  A  DECIMAL  -  DOES  IT  FOLLOW  A  DIGIT 

GO  TO  (37,38),  ITl 

IT  1=2 

LDGIT=1 

IT2=2 

GO  TO  27 

ROUTINE  TO  CONVERT  ALPHABETIC  ARRAY  TO  FLOATING  POINT  NUMBER 
DATA  (NDATA)=GIT(AUXA, 1 , LDGIT, I . )+GIT(AUXB, 1,10,0.) 

DATA  (NDATA)=DATA(NDATA)*SIGN 

IS  ALL  DATA  FURNISHED  YES-RETURN  NO  INCREASE  N  DATA  KEEP  ON 

IF  (NDATA-MAXNO)  41,39,39 

ICC=0 

RETURN 

NDATA=NDATA+1 
GO  TO  25 

FORMAT  (2A1,9A2,60A1) 

FORMAT  (IHl) 

FORMAT  (5X,2A1 ,9A2,60A1) 

FORMAT  (2A1,9A2,60A1) 

FORMAT  (10X,20HC0MMAND  NOT  IN  TABLE) 

END 


FUNCTION  GIT  (TCARD,J,JLAST, SHIFT) 

C  Converts  alphabetic  array  to  floating  point  numbet 


DIMENSION  TCARD(IO),  A(10) 

DATA  A(1)/1H1/,A(2)/1H2/,A(3)/1H3/,A(4)/1H4/,A(5)/1H5/,A(6)/1H6/ 
DATA  A(7)/1H7/,A(8)/1H8/,A(9)/1H9/,A( 10)/1H0/ 

GIT-O. 

TEN=10. 
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DO  3  JNOW=J,JLAST 
TTEST=TCARD(JNOW) 

CHECK  FOR  LAST  ENTRY 
IF  (TTEST.EQ.O.)  GO  TO  4 
FIND  NUMBER  AND  COMPUTE  VALUE 
DO  2  NUMB=1,10 
IF  (TTEST-A(NUMB))  2,1,2 
ZTEST=NUMB 

IF  (ZTEST.EQ.IO.)  ZTEST=0. 

SUM=SUM*TEN+ZTEST 

GO  TO  3 

CONTINUE 

CONTINUE 

IF  (SHIFT)  6,5,6 

FI»JNOW-l 

SUM=SUM*(0.1**FI) 

GIT=SUM 

RETURN 

END 


SUBROUTINE  STHYD 

THIS  SUBROUTINE  STORES  THE  COORDINATES  OF  HYDROGRAPHS. 

COMMON /BLOCKl/  OCFS(30O , 6) ,DATA( 310) ,CFS(300) ,CTBLE( 50 , 1 1 ) , 
&RAIN(300),ROIN(6), 

&A(20,6),Q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 
&ZALFA(20),IEND(6),DA(6),DIST(6),SEGN(6),DT(6),PEAK(6),ISG(6), 
&NPU , NHD , NER , MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 

DIMENSION  DUMMY(300) 

• 

ID=DATA(l) 

NHD=DATA(2) 

DT(ID)=DATA(3) 

IF ( KCODE. EQ.OGO  TO  10 
DATA ( 4 ) =DATA ( 4 ) / 2 . 590 
DO  11  J=5,305 
DATA(J)=DATA(J)/. 02832 
CONTINUE 
DA(ID)=DATA(4) 

J=5 

REMAINING  DATA  ARE  FLOW  RATES 
OCFSd  ,ID)=-DATA(J) 

PEAK(ID)  =  1. 

RO  =  DATA(J) 

DO  4  1=2,300 
J=J+1 

OCFS(I,ID)=DATA(J) 

RO  =  RO  +  OCFSd,  ID) 

IS  FLOW  RECEDING 

IF  (OCFSd, ID)-OCFS(I-l, ID))  1,2,2 
HAS  FLOW  RECEDED  TO  CUTOFF  RATE 
IF  (OCFSd, ID))  5,5,4 
DETERMINE  PEAK  FLOW 


IF(OCFS(I,ID)  -  PEAK(ID))  4,4,3 
PEAK(ID)  -  OCFS(I,ID) 

CONTINUE 

IEND(ID)=I-1 

M-IEND(ID) 

ROIN(ID)  =  (RO*DT(ID))/(DA(ID)*645.333) 

IF(NPU.LE.O)GO  TO  7 
IF(ICODE.EQ.O)GO  TO  6 
ROINl=-ROIN(ID)*25.4 
DA1-DA(ID)*2.590 
PEAK1=PEAK( ID)*. 02832 
DO  13  J=«1,M 

DUMMY( J )=OCFS ( J , ID )*0. 02832 
CONTINUE 

WRITE(7, 14)ID,NHD,DT(ID) ,DA1 ,PEAK1 ,R0IN1 , lEND(ID) , ICODE 
WRITE(7 , 15)(DUMMY(I) ,1=1 ,M) 

RETURN 
PUNCH  CODE 

WRITE(7,8)ID,NHD,DT(ID) ,DA(ID) ,PEAK(ID) ,ROIN(ID) ,IEND(ID) , ICODE 
WRITE  (7,9)  (0CFS(J,ID),J=1,M) 

RETURN 

FORMAT(  'RECALL  HYD' ,T21 ,'ID=' ,11 ,T29,'HYD  N0=' , 13 ,T42 , 'DT=' ,F9 
&6,'  HRS',T61,'DA=',F8.3,'  SQ  MI'/T21 , 'PEAK=' ,F7. 0, 'CFS' ,T40, 'R0=' 
&F6.3,"  INCHES  ",T59,"NO  PTS  =",I3/T21 ,"C0DE=",I1/T21 , 

&"FLOW  RATES") 

FORMAT  (T21,7F8,0) 

FORMAT ("RECALL  HYD" ,T21 ,"ID=" , II ,T29,"HYD  NO  =",I3,T42, 

S."DT="  ,F9 . 6 ,  "HRS"  , T6 1 ,  "DA="  ,F8 . 3 ,  "SQ  KM"/T2 1 ,  "PEAK" ,  F7 . 2 , 
&"CMS",T40,"R0=",F6.0,"  MM  ",T59,"NO  PTS=",I3/T21 ,"CODE=", 
&I1/T21,"FL0W  RATES") 

FORMAT  (T21,7F8.2) 

END 


SUBROUTINE  RECHD 

THIS  SUBROUTINE  RECALLS  PREVIOUSLY  COMPUTED  AND  PUNCHED 
HYDROGRAPHS 

COMMON/BLOCKl/  OCFS(300,6) ,DATA(310) ,CFS(300) ,CTBLE(50 , 1 1 ) , 
&RAIN(300),ROIN(6), 

&A(20,6),Q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 
&ZALFA(20),IEND(6) ,DA(6) ,DIST(6),SEGN(6),DT(6),PEAK(6),ISG(6) , 
&NPU , NHD , NER , MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 

MET1=DATA(8) 

IF(MET1.EQ.0)G0  TO  2 

DATA(4)=DATA(4)/2.590 

DATA(5)=DATA(5)/.02832 

DATA(6)=DATA(6)/25.4 

M=DATA(7) 

DO  3  I=9,M+9 

DATA( I )=DATA( I ) /O . 02832 

CONTINUE 

ID=DATA( 1 ) 

NHD*DATA(2) 


DT(ID)=DATA(3) 

DA(ID)=DATA(4) 

PEAK(ID)=DATA(5) 

R0IN(ID)=DATA(6) 

IENDUD)=DATA(7) 

M=IEND(ID) 

J  =  9 

C  REMAINING  DATA  ARE  FLOW  RATES 

,DO  I  1=1, M 
OCFS(I,ID)=DATA(J) 

1  J=J+1 

RETURN 
END 


SUBROUTINE  CMPHYD 

C  This  subroutine  develops  a  unit  hydrograph,  converts  rainfall  data 
C  into  runoff  by  calling  the  soil  moisture  finite  difference  model, 

C  and  sums  these  two  to  produce  the  storm  runoff  hydrograph. 

COMMON/BLOCKl /  OCFS ( 300 , 6 ) , DATA( 3 1 0 ) , CFS ( 300 ) , CTBLE ( 50 , 1 1) , 
&RAIN(300),ROIN(6), 

&A(20,6),q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 
&ZALFA(20),IEND(6),DA(6),DIST(6),SEGN(6),DT(6),PEAK(6),ISG(6), 
&NPU , NHD , NER , MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 


DIMENSION  DUMMY (300) 
TEMP=0. 


C  Input  data  read  into  subroutine 


ID=DATA( 1 ) 

NHD=DATA(2) 

DT(ID)=DATA(3) 

IF(KCODE.NE.0)THEN 

C  Convert  metric  to  imperial 

DATA ( 4 ) =DATA ( 4 ) / 2 . 590 
IF(DATA(6).LT.0)GO  TO  40 
DATA(  6  )=DATA(  6 )  /  0 . 3048 
DATA(7)=DATA(7)/1.6 

ENDIF 

40  DA(ID)=DATA(4) 

C 


C  Data  items  6  and  7  normally  hold  watershed  height  and  length  and 
C  from  these  the  constants  5ac(recession  constant)  and  Tp(time  to  peak) 
C  can  be  calculated. 

C  If  XK  and  Tp  are  known  however,  they  can  be  entered  instead 


n  n 


s 


C  and  a  negative  sign  is  put  before  their  values. 

IF  (DATA(6).LT.0.)THEN 
XK=-DATA(6) 

TP— DATA(7) 

ELSE 

HT-DATA(6) 

XL=DATA(7) 

SLOPE=HT/XL 

XLDW=(XL**2.)/DA(ID) 

XK=27.0*(DA(ID)**.231)*(SLOPE**(-.777))*(XLDW**.124) 

TP=4.63*(DA(ID)**.422)*(SLOPE**(-.46))*(XLDW**.133) 

ENDIF 


C  The  storm  runoff  array  is  intialised  to  0,  and  peak  of  hydrograph  to  1 

DO  4  1=1,300 

4  0CFS(I,ID)=0. 

PEAK(ID)=1. 

C  Compute  'N'  by  iteration 
XN=5.0 
XKTP=XK/TP 
DO  6  1=1,50 

TINF=1.+SQRT(1./(XN-1.)) 

XN1=.05/(XKTP*(ALOG(TINF/(TINF+.05))+.05))+1. 

DIFF=ABS(XN1-XN) 

IF  (DIFF-.OOl)  7,7,5 

5  XN=XN1 

6  CONTINUE 
WRITE  (6,29) 

29  F0RMAT('  N  DID  NOT  CONVERGE  AFTER  50  ITERATIONS.') 

GO  TO  28 

C  Compute  'Cl' 

7  DELT=TINF/100. 

TC1=0. 

XN1P=XN-1. 

XN1M=1.-XN 

DO  8  1=2,101 
TCl-TCl+DELT 

8  CFS(I)=(TC1**XN1P)*EXP(XN1M*(TC1-1.)) 

SUM=CFS(101  )/2. 

DO  9  1=2,100 

9  SUM=SUM+CFS(I) 

C1=SUM*DELT 


Compute  'B' 

CFSII-CFS(lOl) 

TTINF=TINF*TP 
TREC1*TTINF+2.*XK 
EEE-EXP ( ( TTINF-TREC 1 ) /XK) 

XK1-3.*XK 

B=645.333/(C1+CFSII*(XKTP*(1.-EEE)+EEE*(XK1/TP))) 


c 

C  Compute  'QP'  and  'CFSI' 

C 

QP=(B*DA(ID))/TP 

CFSI=QP*CFS(101) 

CFSR1=CFSI*EEE 
IF(ICODE.EQ.O)GO  TO  45 
QPl-QP*. 02832 
WRITE(6,38)XN,QP1 

38  FORMAT('  Shape  constant,  N  =  ',F6.3/'  Unit  peak  =  ',F10.1,1X 
&,'cms'/) 

GO  TO  44 

45  WRITE  (6,30)  XN,QP 

30  FORMAT('  Shape  constant,  N  =  ',F6.3/'  Unit  peak  *  ',F10.1,1X 
*,'cms'/) 

C 

44  CONTINUE 

C 

C  Determine  the  incremental  runoff 
C 

IF(KCODE.NE.O)THEN 

IF(DATA(8).LT.O)GO  TO  13 

C  Convert  rainfall  data  from  mm  to  inches. 

DO  34  K=8,308 

DATA(K)=DATA(K)/25.4 

34  CONTINUE 
ENDIF 

C 

35  J«8 

IF  (DATA(J))  13,10,10- 

10  RAIN(1)=DATA(J) 

DO  11  1=2,300 

J=J+1 

RAIN(I)=DATA(J.) 

IF  (RAIN(I)-RAIN(I-l))  12,11,11 

11  CONTINUE 

12  NUMB=I-1 

13  CONTINUE 

DO  5555  1=1,300 
5555  DATA(I)=0. 

C 

C 

TEMP=DT(ID) 

C 

CALL  SOILM(TEMP, NUMB, RAIN, DATA) 

C  Subroutine  returns  a  vector  of  runoff  values  from  the  soil  moisture  model 
C  If  no  runoff  has  been  generated  by  the  soil  water  model,  then  the  simulation 
C  stops. 

DO  100  1=1, NUMB 
IF(DATA(I).EQ.0.)G0T0  100 
GOTO  200 
100  CONTINUE 

WRITE (6, 300) 

300  F0RMAT('  Soil  water  model  generated  no  runoff'/ 

&'  Simulation  terminates') 


200  CONTINUE 


Compute  unit  hydrograph 

T2=0. 

CFS(1)=0. 

DO  20  1=2,300 
T2=T2+DT(ID) 

IF  (T2-TTINF)  16,16,17 

16  CFS(I)=QP*((T2/TP)**XN1P)*EXP(XN1M*(T2/TP-1.)) 

GO  TO  20 

17  IF  (T2-TREC1)  18,18,19 

18  CFS(I)=CFSI*EXP((TTINF-T2)/XK) 

GO  TO  20 

19  CFS(I)=CFSR1*EXP((TREC1-T2)/XK1) 

IF  (CFS(I)-l.)  21,21,20 

20  CONTINUE 
1=300 

21  ICND=I 
C 

C 

C  Compute  the  storm  runoff  hydrograph  by  summing  the  unit  hydrograph  and 
C  the  runoff  from  the  soil  moisture  model. 

C 

C 

DO  24  J=2,NUMB 

N=J+ICND-2 

IF  (N-300)  23,23,22 

22  N=300 

23  1  =  2 

DO  24  K=  J,N 

0CFS(K,ID)=0CFS(K,ID)+DATA(J)*CFS(I) 

I=I+l 

24  CONTINUE 
C 

C  Compute  the  runoff  volume  and  determine  the  peak. 

C 

C 

RO  =  0. 

DO  26  I  =  2,N 

RO  =  RO  +  0CFS(I,ID) 

IF  (OCFS(I,ID)-PEAK(ID))26,26,25 

25  PEAK(ID)=  0CFS(I,ID) 

26  CONTINUE 
lEND  (ID)  =  N 

ROIN(ID)=(RO*DT(ID)-)/(DA(ID)  *  645.333) 

C 

C  PUNCH  CODE 

IF  (NPU)  28,28,27 

27  IF(ICODE.EQ.O)GO  TO  39 

ROINl=ROIN(ID)*25.4 
DA1=DA(ID)*2.590 
PEAK1=PEAK(ID)*. 02832 

DO  41  J-1,N 

DUMMY(J)=OCFS(I,ID)*0.02832 
41  CONTINUE 


WRITE(7 , 37 ) ID , NHD , DT( ID) ,DA1 , PEAKl , ROINl , IEND( ID) , ICODE 
WRITE(7,42)(DUMMY(I) ,1=1 ,N) 

RETURN 

39  WRITE(7,31)ID,NHD,DT(ID),DA(ID),PEAK(ID),ROIN(ID),IEND(ID), ICODE 
WRITE  (7,32)  (OCFS(I,ID),I=l,N) 

28  RETURN 

C 

31  FORMAT(  'RECALL  HYD' ,T21 ,'ID=' ,11 ,T29 ,'HYD  NO=' , 13 ,T42 , 'DT=' ,F9. 

&6,'  HRS',T61,'DA=',F8.3,'  SQ  MI'/T21 ,'PEAK=' ,F7 .0 , 'CFS' ,T40 , 'RO=' , 
&F6.3,'  INCHES', T59, 'NO  PTS=' ,I3/T21 ,"CODE=" , II /T21 , 'FLOW  RATES') 

37  FORMATC  'RECALL  HYD' ,T21 ,'ID=' , II ,T29 , 'HYD  NO=' , 13 ,T42 , 'DT=' ,F9. 
&6,'  HRS',T61,'DA=',F8.3,'  SQ  KM'/T21 ,'PEA1C=' ,F7.2,'CMS' ,T40,'RO=' , 
&F6.0,'  MM  ',T59,'NO  PTS=' , I3/T21 ,"CODE=" , II /T21 , 'FLOW  RATES') 

42  FORMAT  (T21,7F8.2) 

32  FORMAT  (T21,7F8.0) 

END 


SUBROUTINE  SOILM( DT , IR , CUMRAIN , DATA) 

C  A  physically  based  parameter  infiltration  model  which  simulates  near  surfa 
C  soil  water  movement,  and  hence  runoff. 

C  Variables  used  in  this  subroutine 


C 

TIME 

Time  when  simulation  begins  (hours). 

c 

SRI 

Soil  water  content  at  saturation  layer 

1. 

c 

SR2 

(m3/m3)  layer 

2. 

c 

SR3 

layer 

3. 

c 

NLA 

Number  of  cells  in  layer  1. 

c 

NLB 

Number  of  cells  in  layer  2. 

c 

NL 

Total  number  of  cells  in  column 

c 

SATCON 

Saturated  permeability  (ms-1)  layer  1. 

c 

SATCON2 

layer  2. 

■ 

c 

SATCON3 

layer  3. 

c 

EMAX 

Maximum  evaporation  during  the  day  (ms- 

1). 

c 

SIMDUR 

Simulation  duration  (hours). 

c 

DETCAP 

Surface  detention  capacity  (m). 

c 

AF 

Simulation  iteration  period  (secs). 

c 

WT 

Write-out  time  period  (hrs). 

c 

THETA 

Initial  soil  water  content  for  each  cell  (m3/m3). 

c 

TCOM 

Thickness  of  each  cell. 

c 

ALR 

Rain  start  time  (hours). 

c 

AMR 

Rain  stop  time. 

c 

NQ 

Number  of  observations  on  suction  moisture  curve. 

c 

X 

Moisture  values ....  layer  1  (m3/m3). 

c 

Y 

Suction  values . layer  1  (bars). 

c 

X2 

layer  2. 

c 

Y2 

layer  2. 

c 

X3 

layer  3. 

c 

Y3 

layer  3. 

c 

IR 

Number  of  rainfall  observations. 

c 

DT 

Rainfall  data  time  increments  (tiours). 

c 

CUMRAIN 

Cumulative  rainfall  data  at  DT  time  increments  (inches). 

c 

NSCOL 

Number  of  soil  columns. 

c 

IPCAREA 

Percent  area  of  soil  column. 
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lOUT 


Determines  amount  of  output. 
1  -  total  output 
0  -  shorter 


Note : 

If  SRI,  SR2,  SR3,  SATCON,  SATC0N2,  SATC0N3,  DETCAP,  THETA,  X,  X2,  or  X3 
are  proceeded  by  an  'A',  then  the  variable  type  Is  double  precision 
rather  than  real.  If  SRI,  SR2,  SR3,  SATCON,  SATC0N2,  SATC0N2,  DETCAP, 
OR  THETA  are  preceeded  by  an  'S', then  the  variable  represents  the 
standard  deviation  of  that  particular  soil  hydrological  characteristic. 

SCURVl  Standard  deviation  of  soil  moisture  curve  for  layer  1 

SCURV2  layer  2 

SCURV3  layer  3 


C  - 

C  INITIAL  SECTION 

C  - 

C 

C 

C 

DIMENSION  FLUX(20),TCOM(20),SWP(20),THETA(20),COND(20) 
DIMENSION  VOL(20) ,ANFLUX(20) ,AVCOND(20) ,DEPTH(20) ,DIST(20) 
DIMENSION  X(20) ,Y(20) ,G(20) ,GZ(20) ,FSWP(20) ,CNT(20) 
DIMENSION  CUMRAIN(25l),Z(20),PPT(250),XP(20),FS(20) 
DIMENSION  DATA(300) ,WDATA(300, 10),HPOT(20) 

DIMENSION  G2(20) ,Y2(20) ,X2(20) ,GZ2(20) ,Z2(20) 

DIMENSION  G3(20) ,Y3(20) ,X3(20) ,GZ3(20) ,Z3(20) 

DIMENSION  RSAT(20) 

DIMENSION  AX(20) ,AX2(20) ,AX3(20),ATHETA(20) 

DIMENS  ION  XNEW( 20 ) , YNEW( 20 ) , X2NEW( 20 ) , Y2NEW( 20 ) , 

&  X3NEW(20) ,Y3NEW(20) 

DOUBLE  PRECISION  G05DDF 
DOUBLE  PRECISION  DLOGIO 

DOUBLE  PRECISION  ATHETA, AX, AX2 , AX3 , ADETCAP , ASRl , ASR2 , ASR3 , 

*  ASATCON , ASATCON2 , ASATC0N3 , BSATCON, BSATC0N2 , BSATC0N3 , 

*  SDETCAP , SSRl , SSR2 , SSR3 , STHETA , SSATCON, SSATC0N2 , SSATCON3 , 

*  SCURVl ,SCURV2,SCURV3 
C 

C 

C 

C  READ  IN  DATA 

C  - 

C 

C 

C 

READ(25, 1000)TIME,ALR,AMR,SIMDUR 
READ(25, 1000)I0UT 
READ(25, 1000)AF,WT 
R£AD(25, 1000)NSC0L 


C  The  array  RAIN  which  is  passed  to  the  subroutine  as  a  cumulative 
C  rainfall  total  is  in  Inches. This  has  to  be  transfered  to  array 
C  PPT  which  is  in  m  and  represents  the  total  for  each  time  increment. 
IRR»IR-1 
DO  100  1=1, IRR 

100  PPT(I)=(CUMRAIN(I+1)-CUMRAIN(I))*.0254 


DO  34543  W=1,NSC0L 

C  For  each  soil  column  in  turn,  read  in  data  and  proceed  through 

C  simulation  to  determine  runoff 

READ( 25 , 1000 ) IPCAREA 
READ(25,1000)NL,NLA,NLB 
READ(25, 1000) (TCOM(l), 1=1, NL) 

READ( 25 , 1000 )EMAX , ADETCAP , SDETCAP 
READ(25,1000)ASR1 ,SSR1,ASR2,SSR2,ASR3,SSR3 

READ( 25 , 1000 ) ASATCON , SSATCON , ASATCON2 , SSATC0N2 , ASATC0N3 , SSATC0N3 
READ(25, 1000 )(ATHETA(I), 1=1, NL) 

READ(25,1000)STHETA 
READ(25, 1000)NQ 
READ(25,1000)(AX(I),I=1,NQ) 

READ(25,1000)(Y(I),I=1,NQ) 

READ(25, 1000)SCURV1 
R£AD(25,1000)(AX2(I),1=1,NQ) 

READ(25,1000)(Y2(I),I=1,NQ) 

READ(25,1000)SCURV2 
R£AD(25, 1000)(AX3(I) ,1=1 ,NQ) 

READ(25,I000)(Y3(I),r=l,NQ) 

R£AD(25, 1000)SCURV3 
1000  FORMAT (V) 

NQJ=NQ 

NLL=NL+1 

1F(AMR.LT.ALR)THEN 

AMR=AMR+24.0 

ENDIF 

C 

C 

C  CHECK  DATA  INPUTS 

C  - 

C 

NERROR=0 

C  Check  number  of  cells  in  soil  column 
IF(NLA+NLB.GE.NL)THEN 
WRITE(6, 1015) 

1015  F0RMAT('  Error-NLA,NLB,NL') 

NERR0R=NERR0R+1 

ENDIF 

C 

C  Check  dimensions  of  input  vectors 

IF(NQ.GT.20.0R.NL.GT.20.0R.IR.GT.250)THEN 
WRITE(6, 1020) 

1020  F0RMAT('  Error-limit  exceeded ,NQ , NL, IR' ) 

NERR0R=NERR0R+1 


C  Check  rainfall  passed  from  CMPHYD 
KN=IR-1 
DO  50  1=1, KN 

IF ( CUMRAIN ( 1+ 1 ) . LT . CUMRAIN ( I ) )THEN 
WRITE (6, 1030) 

1030  FORMAT('  Error-not  cumulative  rainfall  totals') 

NERR0R=NERR0R+1 

ENDIF 

50  CONTINUE 
C 

C  Check  that  initial  moisture  content  of  each  cell  lies  within  the  range  of 
C  the  suction  moisture  curve  and  does  not  exceed  stated  saturated  moisture 
C  content. 

DO  51  1=1, NLA 

IF(ATHETA(I).GT.ASR1)THEN 

WRITE(6,1050) 

1050  F0RMAT('  Error-THETA  larger  then  sat  moisture  concent(l)') 

NERR0R=NERR0R+1 
ENDIF 

I F  ( ATHE TA  < I ) . GT . AX ( NQ ) . OR . ATHETA ( I ) . LT . AX ( 1 ) ) THE  N 
WRITE(6,1055) 

1055  F0RMAT('  Error-THETA  outside  range  of  curves-(l)') 

ENDIF 

51  CONTINUE 
NLAA=NLA+1 
NLH=NLA+NLB 

DO  52  I=NLAA,NLH 

IF(ATHETA(I).GT.ASR2)THEN 

WRITE(6,1060) 

1060  F0RMAT('  Error-THETA  larger  than  sat  moisture  content(2)') 

NERR0R=NERR0R+1 
ENDIF 

I F ( ATHET A ( I ) . GT . AX2 ( NQ ) . OR , ATHETA( I ) . LT . AX2 ( 1 ) )THEN 
WRITE(6, 1065) 

1065  F0RMAT('  Error-THETA  outside  range  of  curve-(2)') 

NERR0R=NERR0R+1 

ENDIF 

52  CONTINUE 
NLBB=NLB+NLA+1 
DO  53  I=NLBB,NL 

IF(ATHETA(I).GT.ASR3)THEN 

WRITE(6,1070) 

1070  F0RMAT('  Error-THETA  larger  than  sat  moisture  content(3)') 

STOP 

ENDIF 

IF( ATHETA( I ) .GT . AX3 ( NQ ) .OR. ATHETA( I ) . LT. AX3 ( 1 ) )THEN 
WRITE(6,1075) 

F0RMAT('  Error-THETA  outside  range  of  curve  -(2)') 


1075 


o  o  n  n  ui 


NERR0R=NERR0R+1 
ENDIF 
3  CONTINUE 

IF  (NERROR.NE.O)THEN 

WRITE (6, 1076 )NERROR 

076  FORMAT('  SOILM:  number  of  input  data  errors  ',12, 

&'Simulation  terminates') 

STOP 

ENDIF 


C 

C  DEPTH  CALCULATION 

C  - 

C 

C' 

C 

C  The  variable  DEPTH  is  calculated.  This  refers  to  the  distance  from 
C  ground  level  to  any  cell  midpoint. 

C  DIST  refers  to  the  distance  between  any  two  adjacent  cell  midpoints. 

C 

DIST(l)=TCOM(l)/2. 

DEPTH(1)=*DIST(1) 

DO  110  1=2, NL 

DEPTH( I )=DEPTH( I-l )+0 . 5*(TCOM( I-l )+TC0M( I ) ) 

110  DIST(I)=0.5*(TCOM(I-1)+TCOM(I)) 

C 

C 

C 

C 

C  PARAMETER  VARIABILITY 

C  - 

C 

C  Five  input  variables,  detention  capacity,  soil  water  content  at 
C  saturation,  soil  moisture  content  at  given  tensions,  saturated  conductivity 
C  and  initial  moisture  content  are  varied  stochastically. 

C  NAG  functions  are  called  which  return  a  'psuedo  random'  value  from  a 
C  distribution  with  a  given  standard  deviation  and  mean. 

C  All  are  assumed  to  have  a  normal  distribution  except  the  saturated 
C  conductivity  which  takes  on  a  lognormal. 

C 

C 

C 

C  Generate  only  one  set  of  stochastic  variables  to  run  in  HYMO. 

C 

C 

C  RANDOM  PARAMETER  VALUE 

r  - 


non  •—  ooi—  on 


C 

WRITE(6,1079) 

1079  F0RMAT('  INCREMENTAL  RUNOFF-Parameter  variability  included'//) 

Detention  capacity. 

DETCAP=G05DDF(ADETCAP,SDETCAP) 

IF(DETCAP.LT.0.)DETCAP=0.0 

SD=SDETCAP 

WRITE(6,1180)SD 

180  F0RMAT('  SD  of  detcap  ',F5.3) 

Soil  water  content  at  saturation 
SR1=-G05DDF(ASR1,SSR1) 

SR2-G05DDF(ASR2,SSR2) 

SR3-G05DDF(ASR3,SSR3) 

SDl-SSRl 

SD2-SSR2 

SD3=SSR3 

WRITE(6, 1181 )SD1,SD2,SD3 

181  F0RMAT('  SD  of  saturated  soil  content' ,F5. 3 , '  layer  1'/ 

&  '  ',F5.3,'  layer  2'/ 

&  '  ',F5.3,'  layer  3') 

Soil  moisture  content  at  given  tensions 
Layer  1 

CALL  SMCURV ( S R 1 , NQ , AX , Y , XNEW , YNE W , SCURV 1 ) 

DO  120  1-1,20 
X(I)=XNEW(I) 

120  Y(I)-YNEW(r) 

C  Layer  2 

CALL  SMCURV( SR2 , NQ , AX2 , Y2 , X2NEW , Y2NEW, SCURV2 ) 

DO  130  1=1 ,20 
X2(I)=X2NEW(I) 

130  Y2(I)=Y2NEW(I) 

C  Layer  3 

CALL  SMCURV ( SR3 , NQ , AX3 , Y3 , X3NEW , Y3NEW , SCURV3 ) 

DO  140  1-1,20 
X3(r)-X3NEW(I) 

140  Y3(I)-Y3NEW(I) 

SDl-SCURVl 

SD2-SCURV2 

SD3-SCURV3 

WRITE(6, 1 182)SD1 ,SD2,SD3 

1182  F0RMAT('  SD  of  suction  moisture  curve',  F5.3,'  layer  1'/ 

4  '  ',  F5.3,'  layer  2'/ 

4  '  ',F5.3,'  layer  3') 

C 

C  Saturated  conductivity  for  each  layer 
BSATCON-DLOG 1 0 ( ASATCON ) 

SATC0N-G05DDF( BSATCON , SSATCON ) 


SATCON=10**SATCON 
BSATC0N2=DL0G 1 0 ( ASATC0N2 ) 

SATC0N2=G05DDF ( BSATC0N2 , SSATC0N2 ) 

SATCON2=10**SATCON2 
BS ATCON  3=DL0G 1 0 ( AS ATCON3 ) 

SATC0N3=G05DDF( BSATCON3 , SSATCON3 ) 

SATC0N3-10**SATC0N3 

SDl=SSATCON 

SD2=*SSATCON2 

SD3=SSATCON3 

WRITE(6,U83)SDl,SD2,SD3 

1183  FORMAT('  SD  of  sat  conductivity' ,F5. 3 , '  layer  1'/ 

&  '  ',F5.3,'  layer  2'/ 

&  '  ',F5.3,'  layer  3') 

C 

C  Initial  moisture  content 
DO  150  1=1, NL 

150  THETA(I)=G05DDF(ATHETA(I),STHETA) 

C  Check  on  Initial  soil  moisture  values 

DO  160  1=1, NLA 

IF(THETA(I).GE.X(20))THETA(I)=X(20)-0.001 
160  IF(THETA(I).LE.X(1))THETA(I)=X(1)+0.001 

DO  170  I=NLAA,NLH 

IF(THETA(I).GE.X2(20))THETA(I)=X2(20)-0.001 
170  IF(THETA(I).LE.X2(1))THETA(I)=X2(1)+0.001 

DO  180  I=NLBB,NL 

IF(THETA(I).GE.X3(20))THETA(I)=X3(20)-0.001 
180  IF ( THETA( I ) . LE . X3 ( 1 ) )THETA( I )=X3 ( 1 )+0 . 00 1 

SD=STHETA 
WRITE(6,118A)SD 

1184  F0RMAT('  SD  of  initial  water  content' ,F5. 3) 

C 

C 

C 

C  HYDRAULIC  CONDUCTIVITY  CALCULATION 

C  - 

C 

c 

c 

C  The  hydraulic  conductivity  is  calculated  from  suction  moisture 
C  data  for  each  layer, 

NQJ=NQ 

CALL  HYDC0N(X,SATC0N,SR1,Z,Y) 

CALL  HYDC0N(X2 , SATC0N2 , SR2 , Z2 , Y2 ) 

CALL  HYDCON ( X3 , SATC0N3 , SR3 , Z3 , Y3 ) 

C 
C 
C 
C 
C 


WRITE-OUT  INITIAL  CONDITIONS 


o  o 


c 

C  Write-out .suction  moisture  curve  and  generated  K-values. 

C 

WRITE (6, 1080) 

1080  FORMATC'OGENERATED  K-MOISTURE  CURVE'/ 

&'  Millington-Quirk  Method'/ 

&'  Layer  1' ,26X, 'Layer  2' ,26X,'Layer  3'/ 

&3('  Moisture  Suction  Unsat  K  ')) 

DO  175  I=*l,20 

175  WRITE(6,1090)X(I),Y(I),Z(I),X2(I),Y2(I),Z2(I),X3(I),Y3(I),Z3(I) 

1090  FORMATdH  ,  3(F6. 3 , 2X,F8. 3 ,F15. 12 , 2X) ) 

Write-out  start  conditions. - 

WRITE(6,1100) 

1100  FORMAT ('OSTART  CONDITIONS  '/) 

WRITE(6,1110)TIME 

1110  FORMAT('  Simulation  start  time' ,F4. 1 ,'hrs') 

WRITE(6, 1130)ALR,AMR 

1130  FORMAT('  Precipitation  begins  at  ' , F4. 1 , 2X, 'and  ends  at  ',F4.1) 
WRITE(6,U^0)DT 

1140  F0RMAT('  Rainfall  data  time  increment  =  ' ,F6.4 , 2X, 'hrs' ) 
WRITE(6,1120)AF 

1120  F0RMAT('  Time  increment  for  iteration  period  =  ',F6.1, 

&2X,'secs'/) 

WRITE (6,1150 )EMAX , DETCAP 

1150  F0RMAT('  Maximum  evaporation  during  the  day  *  ' ,F10.8,2X,'ms-l'/ 

&'  Surface  detention  capacity  =  ',F6.4,2X,'m'//) 

C 

C  Calculate  initial  relative  saturation  of  each  cell  in  soil  column 
DO  1151  1=1, NL 

IF(I.LE.NLA)RSAT(I)=THETA(I)/SR1 

IF(I.GT.NLA.AND.I.LT.NLBB)RSAT(I)=THETA(I)/SR2 

IF(I.GE.NLBB)RSAT(I)=THETA(I)/SR3 

1151  CONTINUE 

WRITE(6,1152) 

1152  F0RMAT('  INITIAL  SOIL  COLUMN  CONDITIONS'//) 

WRITE(6,1153) 

1153  F0RMAT(11X,'SAT',8X,'SAT  HYD' ,6X, 'CELL' , IX, 'DEPTH' , 

&2X , ' INITAL' , 2X , ' REL' / 

&1H  , lOX, 'THETA' ,7X,'COND',9X, 'NO', lOX, 'THETA' ,2X, 'SAT'/ 

6ilH  , 10X,'m3/m3',7X,'ms-l',14X,'m' ,5X,'m3/m3'/) 

WRITE (6,1154)SR1, SATCON , DE  PTH( 1 ) , THETA( 1 ) , RS AT ( 1 ) 

1154  F0RMAT('  Layer  1  ' ,F7. 4, 1X,F15. 12 , 3X, ' 1 ' , 2X, F6. 4, IX, F7. 4 , IX, F5. 3) 

IF(NLA.GT.  DTHEN 

DO  1155  1-2, NLA 


n  n  o  o  o  o  n 


WRITE ( 6 , 1 1 56 ) I , DE  PTH ( 1 ) , THETA (I ) , RSAT ( I ) 

1156  FORMATdH  ,  34X,I2 , 2X,F6.4 ,  IX, F7 .4 ,  IX, F5. 3) 

1155  CONTINUE 

ENDIF 

WRITE(6,1157)SR2,SATCON2,NLAA,DEPTH(NLAA),THETA(NLAA),RSAT(NLAA) 

1157  FORMAT('  Layer  2  ' ,F7.4, 1X,F15. 12,2X, 12 , 2X,F6.4 , 1X,F7.4 , 1X,F5. 3) 

IF(NLB.GT.1)THEN 

DO  1158  I=NLA+2,NLH 

WRITE ( 6 , 1 1 59 ) I , DEPTH( I) , THETA( I ) , RSAT( I ) 

1159  FORMATdH  ,34X,I2,2X,F6.4, 1X,F7.4, 1X,F5.3) 

1158  CONTINUE 
ENDIF 

WRITE(6,1160)SR3,SATCON3,NLH+1,DEPTH(NLH+1) ,THETA(NLH+1) , 
&RSAT(NLH+1) 

1160  FORMAT('  Layer  3  ' ,F7 . 4 , 1X,F15. 12 , 2X, 12 , 2X, F6. 4 , IX, F7 . 4 , 1X,F5. 3) 

IF((NL-NLH).GT.1)THEN 

DO  1161  I=NLH+2,NL 

WRITE(6, 1162)I,DEPTHd),THETAd),RSATd) 

1162  FORMATdH  ,  34X,  12 , 2X,F6.4 , 1X,F7 . 4 ,  IX, F5 . 3 ) 

1161  CONTINUE 
ENDIF 

C 


INITIALISATION  OF  VARIABLES 


DO  184  1=1,300 

184  WDATAd,W)=»0.0 
WATI=0.0 

MMM=2 

DO  185  1=2, NL 

185  ANFLUXd)=0.0 
CTIME=TIME*3600 
SRAINI=0.0 
CUMDRN=0. 
CINFIL=0. 

SUMD=0. 

ICOUNT  =0 

BR=AMR-ALR 

EVAPI=0.0 

SOG=THETAd)/SRl 

RTOT=0.0 

ANFILT=0.0 

PPTT=0.0 

TG=0.0 


o  n  o  n  o  o 


c 

c 

c 

c 

c 

c 

C  A  ca 
C  The 
C 


BALANCE  CHECK 


ilculatlon  for  the  water  balance  check, 
initial  soil  water  content  of  the  soil  column. 

DO  190  1=1, NL 

WATI=TC0M( I ) *THETA( I )+WATI 


CURVE  GRADIENTS 


C  Calculations  of  the  gradients  of  the  suction-moisture  curve  and  the 
C  K-moisture  curve  for  each  layer. 

C 

CALL  GRAD(G,GZ,Y,X,Z) 

CALL  GRAD(G2,GZ2,Y2,X2,Z2) 

CALL  GRAD(G3,GZ3,Y3,X3,Z3) 


DYNAMIC  SECTION  -  SIMULATION 


C  This  loop  is  completed  for  each  time  increment  until  end  of  simulation 
C 

ITMAX=SIMDUR*3600/AF 
DO  9995  II=1,ITMAX 
IC0UNT=IC0UNT+AF 
TG=TG+AF 
T=II 


CALCULATE  WATER  VOLUME  OF  EACH  CELL 


DO  200  1=1 ,NL 

V0L( I)=TC0M( I )*THETA( I ) 


V 


'•  **'%**% 
-  •  K  *- 

■ 
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24-hour  clock 


C  Calculate  REAL  TIME  for  current  Iteration  period  using  the  24-hour  clock 
C 

CTIME=CTIME+AF 
IF  (CTIME.GE. 86400 )THEN 
CTIME=CTIME-86400 
ENDIF 
C 
C 

SWP,HP0T,C0ND  CALCULATIONS 


rV--.  AV. 


m 


'**•«•*  X  - 
'  «  ■  k  ' 

itLi 


Calculate  the  soil  water  pressure,  hydraulic  potential  and  conductivity 
for  each  cell  as  conditions  change  during  the  simulation. 

CALL  TW0( 1 , NLA , THETA , X , SWP , Y , G , HPOT , DEPTH , GZ , COND , Z ) 

CALL  TWO ( NLAA , NLH , THETA , X2 , SWP , Y2 , G2 , HPOT , DEPTH , GZ 2 , COND , Z2 ) 

CALL  TWO ( NLBB , NL , THETA , X3 , SWP , Y3 , G3 , HPOT , DEPTH , GZ 3 , COND , Z3 ) 


.V. 


DETERMINE  RAINFALL 


K.m 


Determine  rainfall  per  second  at  end  of  the  current  iteration 
period . 

T1  is  the  time  in  hours  when  the  current  iteration  period  ends. 

Check  that  Tl  is  between  the  rain  start  and  stop. 

If  it  is,  decide  which  element  of  PPT  array  the  data  is  to  be  taken  from 
and  make  SRAIN  equal  to  that  precipitation  per  second. 

If  it  is  not  within  the  storm  period,  set  SRAIN  to  0. 


T1=T*AF/ 3600.0 

IF(T1. LE.(ALR-TIME). OR. T1.GT,( AMR-TIME) )THEN 
SRAIN=0.0 
ELSE 

T2=T1-(AF/3600.) 

IELEM=((T2-(ALR-TIME))/DT)+1 

SRAIN=PPT(IELEM)/(DT*3600.0) 

ENDIF 


c 

c 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 


Increment  precipitation  total  by  amount  of  precipitation  in  current 
iteration  period. 

PPTT*PPTT+( SRAIN*AF ) 


AVERAGE  HYDRAULIC  CONDUCTIVITY 


Average  hyraulic  conductivity  “"o  flow  through  boundary  between 
adjoining  cells  is  weighted  act  .ding  to  its  thickness. 


DO  210  1=2, NL 

210  AVC0ND(I)=(C0ND(I-1)*TC0M(I-1)+C0ND(I)*TC0M(I)) 

&/(TC0M(I-l )+TC0M(I)) 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


BOTTOM  BOUNDARY  CONDITION 


Determine  the  bottom  boundary  condition  under 
water  is  flowing  out  of  the  soil  column  under 


the  assumption  that 
gravity. 


FLUX(NLL)=COND(NL) 


C 

C 

c 

c  FLUX  BETWEEN  CELLS 

- 

C 

C 

C 

C  The  flux  between  each  cell  then  follows  Darcy's  law  in  discrete  form 
C 


220 

C 

C 

c 

c 


DO  220  1=2, NL 

FLUX(I)=(HPOT(I-l )-HP0T(I))*AVC0ND(I)/DIST(I) 


DETERMINE  TOP  BOUNDARY  CONDITIONS 


ono  noo  ooonnon 
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Calculate  the  infiltration  capacity. 

BNCAP=(0 . 0-HPOT( 1 ) )*0 . 5*( SATC0N+C0ND( 1 ) ) /DIST( 1 ) 

Calculate  precipitation  excess 

IF( SPAIN  1 . EQ . SPAIN )THEN 

SUMD= ( SPAIN-ANF I LT ) *AF+SUMD 
ELSE 

SUMD=0.0+SUMD 
ENDIF 

SPAIN1=SRAIN 

Calculate  amount  detained  on  the  surface. 

IF(SUMD.LT.O.O)THEN 
DETAIN=0.0 
ELSE 

DETAIN=SUMD 
ENDIF 

Calculate  evaporation,  the  flux  into  cell  1  and  runoff. 

IF(SPAIN.GT.O.O)  THEN 
C 

EVAP=  0.0 
C 

IF( SPAIN. LT.BNCAP. AND. DETAIN. LE.0.0)THEN 
ANFILT=SRAIN 
ELSE 

ANFILT=BNCAP 
ENDIF  . 

FLL'X(  1  )=ANFILT 
C 

IF(DETAIN.GT.DETCAP)THEN 

SUMD=DETCAP 

DETAIN=DETCAP 

PUNOFF=0.0 

IF(SPAIN.GT.BNCAP)RUNOFF=(SPAIN-BNCAP)*AF 

RTOT=RTOT+RUNOFF 

ELSE 

RUNOFF =0.0 
ENDIF 


C 


RUNOFF =0.0 


n  o 


IF(CTIME.GT. 64300. AND. CTIME.LE.21600)THEN 
EVAP-EMAX/100. 

ELSE 

EVAP=EMAX*SIN ( 2 . *3 . 1 4 1 59*( CTIME-2 1600.) /86400 . ) 

ENDIF 
C 

IF(DETAIN.LE.O.)THEN 
ANFILT=0.0 
FLUX(l)=EVAP*(-l.) 

ELSE 

ANFILT=BNCAP 
FLUX(1)=ANFILT 

DETAIN=DETAIN-(EVAP*AF) 

ENDIF 
C 

ENDIF 
C 
C 

C  CHANGES  IN  SOIL  MOISTURE  CONTENT 

C  - 

C 
C 
C 

SWP(NLL)=-102.0 
DO  230  1=1, NL 

If  SWP  in  cell  is  greater  then  0,  it  is  saturated  and  flux  must 
therefore  be  0. 

IF(SWP(I+1).GE.0.0)FLUX(I+1)=0.0 
C  ANFLUX  represents  the  net  change  in  moisture  content  in  the  cell. 
ANFLUXC I )=FLUX( I )-FLUX( I+l ) 

ANFLUX( I )=ANFLUX( I )*AF 

C  Recalculate  theta  according  to  the  change  influx(per  unit  area). 
THETA( I )=( V0L( I )+ANFLUX( I ) ) /TCOM( I ) 

C  Due  to  recalculation,  theta  may  be  greater  than  possible  water  content 

C  at  saturation  and  therefore  it  is  necessary  to  reset  SWP  to 

C  0  and  theta  to  the  water  content  at  saturation,  the  value  of  which  is 
C  entered  into  the  model. 

IF  (THETA(I).GE.SR1.AND.I.LE.NLA)SWP(I)=0.0 
IF  (THETA(I).GE.SR2.AND.I.GT.NLA.AND.I.LE.NLH)SWP(I)=0.0 
IF(THETA(I).GE.SR3.AND.I.GT.NLH)SWP(I)=0.0 
IF(THETA(I).GE.SR1.AND.I.LE.NLA)THETA(I)=SR1 
IF(THETA(I).GE.SR2.AND.I.GT.NLA.AND.I.LE.NLH)THETA(I)=SR2 
230  IF(THETA(I).GE.SR3.AND.I.GT.NLB)THETA(I)=SR3 
C 
C 
C 
C 


CALCULATE  CUMULATIVE  TOTALS 
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CUMDRN=CUMDRN+FLUX( NLL)*AF 
EVAPI=EVAP*AF+EVAPI 
CINFIL=CINFIL+ANFILT*AF 
S0G=THETA(1)/SR1 


TERMINAL  SECTION  WRITE  OUT 


To  print  out  data  for  every  time  increment  for  which  PPT  data  is 
entered,  check  ICOUNT  to  see  if  that  period  has  passed  by. 
IF(ICOUNT.LT.(DT*3600))  GOTO  9995 
IC0UNT=0 


CALCULATE  TIME  FROM  THE  START 


T=T*AF/3600 

WRITE(6,1170)T 

1170  FORMAT('OSOIL  COLUMN  CONDITIONS  ' ,F7.3, lX.'HRS  SINCE  RAIN  BEGAN'/) 
IFCTG.EQ. 86400. 0)TG=0.0 
C 
C 

C  WRITE-OUT  CONDITIONS  OF  SOIL  COLUMN 

C  - 

C 

C 

IF(IOUT.Eq.O)GOTO  305 
WRITE(6,7780) 

7780  FORMAT('  Cell  Depth  SWP  Theta  Hyd  cond  Net', IX, 
&'flux  Rel  sat') 

DO  300  1=1, NL 

IF(I.LE.NLA)S0G=THETA(I)/SR1 

IF ( I . GT . NLA. AND . I . LT . NLBB ) S0G=THETA( I ) /SR2 

IF(I.GE.NLBB)SOG=THETA<I)/SR3 


300  WRITE ( 6 , 1 1 90 ) I , DE  PTH ( I ) , S WP ( I ) , THETA ( I ) , COND ( I ) , ANFLUX ( I) , SOG 

1190  F0RMAT(I6,3F8.4,2F14.9,F9.3) 

C 

C 

C  WATER  BALANCE  CHECK 

Q  - 

C 

C  Philips  (1964)  simple  water  balance; 

C - 

C 

C 

C  Amount  added 

C  (Initial  soll)-(Current  soil)  =  by  -  Evaporation-  Drainage 

C  (  moisture  )  (  moisture  )  infiltration  loss  loss 

C 

305  WATN=0. 

DO  310  1=1, NL 

310  WATN=TC0M(I)*THETA(I)+WATN 

BAL=WATN-WAT I -C INF I L+E VAP I +CUMDRN 
WRITE ( 6, 1200) BAL 

1200  F0RMAT('0Balance  check  on  soil  column  water  status  =',F12.7) 
BAL=(BAL*100.)/WATN 
WRITE(6,1210)BAL 

1210  FORMAT('  Balance  check  as  column  water  vol.  =',F12.7,'  X'/) 

C 

C 

IF(IOUT.EQ.O)GOTO  306 

WRITE(6, 1220)EVAPI,PPTT,CINFIL,CUMDRN 

1220  F0RMAT('  Cumulative  evaporation  =  ',F12.8/ 

&'  Cumulative  precipitation  =  ',F8.4/ 

Cumulative  infiltration  =  ',F10.6/ 

&'  Cumulative  drainage  =  ',F10.6A) 

306  IF(DETAIN.EQ.DETCAP)THEN 

WRITE(6,1222) 

1222  F0RMAT('  Detention  capacity  exceeded') 

WRITE(6, 1230)RTOT,RTOT/.0254,T 

1230  F0RMAT('  Runoff  total  in  the  last  period' , F 10 . 7 , 2X, 'm' / 

&  '  Runoff  total  in  the  last  period' , FIO. 7 , 2X, 'ins' , 

$  F7.3/) 

ELSE 

WRITE(6,1221)DETAIN 

1221  F0RMAT('  Surface  water  =  ',F10.6) 

WRITE(6,1226) 

F0RMAT('  No  runoff') 


1226 
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C  CREATION  OF  ARRAY  DATA 

C  - 

C 

C 

C  Runoff  is  recorded  in  array  WDATA 

C  The  runoff  for  each  soil  column  is  weighted  according  to  the 
C  percentage  area  which  it  occupies  in  the  catchment  area 
C 

WDATA(MMM,W)=(RTOT/.0254)*(IPCAREA/100. ) 

RTOT=0.0 
MMM=MMM+1 
9995  CONTINUE 


C  End  of  simulation  of  single  soil  column,  it  more  than  one,  then  return  to 
C  to  the  beginning  of  this  subroutine  to  repeat  for  next  soil  column 

34543  CONTINUE 


DO  76567  1=1, MMM 

C  Sum  the  weighted  runoff  for  each  soil  column  to  derive  total  runoff 
C  passed  back  to  CMPHYD  as  DATA 
CUMDATA=0. 

DO  54345  J=1,NSC0L 

CUMDATA=WDATA( I , J)+CUMDATA 
54345  CONTINUE 

DATA(I)=CUMDATA 
76567  CONTINUE 


IR=MMM-l 

RETURN 

END 


SUBROUTINE  HYDC0N(X,SATC0N,SR,Z,Y) 

C  This  subroutine  calculates  hydraulic  conductivity  for  each  layer 
C  from  the  given  soil  moisture  characteristic  curve. 

C  Uses  the  Millington  and  Quirk  method 

DIMENSION  X(20) ,Y(20) ,Z(20) 

DO  845  1=1,20 

IIJ=20-I+l 

XII=X(HJ) 

TOPS»0. 

BOTS=0. 


o  o 


JF-20-J+1 

YJJ=Y(JF) 

8A6  BOTS=((2*J-l)*YJJ**(-2))+BOTS 

II=I 

DO  847  J=II,20 

JF=20-J+1 

YJJ=Y(JF) 

847  TOPS=((2*J+l-2*I)*YJJ**(-2))+TOPS 

JT=20-I+1 

845  Z(JT)=SATCON*(X(II)/SR)*TOPS/BOTS 
RETURN 
END 


SUBROUTINE  TWO ( NA , NB , THETA , X , SWP , Y , G , HPOT , DEPTH , GZ , COND , Z ) 


C  This  subroutine  calculates  soil  water  pressure,  hydraulic  potential 
C  and  hydraulic  conductivity  for  each  cell  as  conditions  change 
C  during  simulation. 


DIMENSION  THETA(20),X(20),SWP(20),Y(20),G(20),HPOT(20), 
&DEPTH(20) ,GZ(20) ,C0ND(20) ,Z(20) 

DO  15  I=NA,NB 
DO  16  J=l,19 

IF(THETA(I).GE.X(J).AND.THETA(I).LT.X(J+1))SWP(I)=Y(J)+G(J)* 

&  (THETA(I)-X(J)) 

16  CONTINUE 

HPOT(I)=SWP(I)-DEPTH(I) 

DO  17  J=l,19 

IF(THETA(I).GT.X(J).AND.THETA(I).LE.X(J+l))COND(I)"Z(J)+GZ(J) 
&  (THETA(I)-X(J)) 

17  CONTINUE 

15  CONTINUE 
RETURN 
END 


SUBROUTINE  GRAD(G,GZ, Y,X,Z) 


C  This  subroutine  calculates  the  gradients  of  the  suction-moisture 
and  hydraulic  conductivity-moisture  curves. 

DIMENSION  G(20) ,GZ( 20 ) , Y(20) ,X(20) ,Z( 20) 

DO  261  1=1,19 

G(I)-(Y(I+1)-Y(I))/(X(I+1)-X(I)) 

261  GZ(I)=(Z(I+1)-Z(I))/(X(I+1)-X(I)) 

RETURN 


1  V"*.” 


SUBROUT INE  SMCURV ( S R , NQ , AX , Y , XNEW , YNEW , SCURV ) 

C  Generates  a  stochastic  suction  moisture  curve  to  be  fed  into 
C  soil  moisture  model 


DOUBLE  PRECISION  G05DDF 
DOUBLE  PRECISION  AX, SCURV 

DIMENSION  AX(20),X(20),XNEW(20),YNEW(20),G(20) ,Y(20) 


C  Determine  the  stochastic  values  of  moisture 
C 

X( 1 )=G05DDF(AX( 1 ), SCURV) 

IF(X(l).LT.O. )X(1)=0.001 
C 

DO  100  1=2, NQ 
X(I)=G05DDF(AX(I) , SCURV) 

100  IF(X(I).LE.X(I-1))X(I)=X(I-1)+0.001 

IF(X(NQ).GE.SR)SR=X(NQ)+0.001 
C 

C  Calculate  gradients  of  this  new  suction-moisture  curve 
c 

NNq=NQ-l 
DO  200  1=1, NNQ 

200  G(I)=(Y(I+1)-Y(I))/(X(I+I)-X(I)) 

C 

C  Calculate  max  and  min  moisture  values,  and  determine  the  size  of 
C  equal  intervals. 

C 

XMAX=RMAX(X,NQ) 

XMIN=RMIN(X,NQ) 

XINT=(XMAX-XMIN)/19. 

C 

C  Determine  the  new  values  of  moisture-equal  intervals 
C 

XNEW(1)=XMIN 
DO  300  1=2,19 

300  XNEW(I)=XNEW(1)+(XINT*(I-1)) 

XNEW(20)=XMAX 

C 

C  Determine  the  associated  new  values  of  suction 
C 

DO  350  1=1,19 
DO  400  J=l,NNq 

IF(XNEW(I).GE.X(J).AND.XNEW(I).LT.X(J+1)) 
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&  YNEW(I)=Y(J)+G(J)*(XNEW(I)-X(J)) 
400  CONTINUE 
350  CONTINUE 

YNEW(20)=Y(NQ) 


RETURN 

END 


FUNCTION  RMAX  (X.NQ) 


Determines  the  maximum  real  in  an  array 


DIMENSION  X(NQ) 

RMAX=X( 1 ) 

DO  10  1=2, NQ 

0  IF(X(I).GT.RMAX)RMAX=X(I) 

RETURN 

END 


FUNCTION  RMIN(X,NQ) 


C  Determines  minimum  real  in  an  array 

DIMENSION  X(NQ) 

C 

RMIN=X( I ) 

DO  10  1=2, NQ 

10  IF(X(I).LT.RMIN)RMIN=X(I) 

C 

RETURN 

END 


SUBROUTINE  PRTHYD 

C  THIS  SUBROUTINE  PRINTS  THE  COORDINATES  OF  A  HYDROGRAPH. 

COMMON/BLOCKl/  OCFS(30O,6) ,DATA(310) ,CFS(300) ,CTBLE( 50 , 1 1 ), 
&RAIN(300),R0IN(6), 

iA(20,6) ,Q(20,6) ,DEEP(20,6) ,ITBLE(50,2),DP(20) ,SCFS(20) ,C(20) , 
4ZALFA(20),IEND(6),DA(6),DIST(6),SEGN(6),DT(6),PEAK(6),ISG(b), 


&NPU , NHD , NER , MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 
DIMENSION  DUMMY (300) 

C  Input  data  is  read  into  the  subroutine. 

ID=DATA( 1 ) 

NPK=DATA(2) 

C  DETERMINE  TYPE  OF  HYDROGRAPH 
IF  (NHD-100)  6,6,2 

1  WRITE  (6,14)  NHD 
GO  TO  7 

2  IF  (NHD-300)  3,3,4 

3  WRITE  (6,15)  NHD 

GO  TO  7 

4  IF  (NHD-500)  1,1,5 

5  WRITE  (6,16)  NHD 

GO  TO  7 

6  WRITE  (6,17)  NHD 

C  POSITIVE  NPK  MEANS  PRINT  ONLY  PEAK  AND  VOLUME 

7  IF  (NPK)  8,8,25 

8  J-0 
M=IEND(ID) 

TIME 1 “TIME 

C  BUILD  TIME  ARRAY  IN  DATA 

DO  9  1=1 ,M 
DATA  (I)=TIME1 

9  TIMEI=TIME1+DT(ID) 

M4=M+4 

M5=M4/5 
DO  22  1=1, M 

DUMMY(I)=OCFS(I,ID)*0.02832 
22  CONTINUE 


WRITE(6,27) 

24  J=J+1 

WRITE ( 6 , 39 ) ( DATA( I ) , DUMMY( I ) , I=J ,M , M5 ) 
IF(J-M5)24,25,25 

25  ROINl=ROIN( [D)*25.4 
PEAK1=PEAK(ID)*0. 02832 
WRITE(6,26)R0IN1 ,PEAK1 

30  IF(NPU)13, 13, 12 

12  WRITE  (7,21)  ID, NPK 

13  RETURN 
C 

14  FORMAT  ( 1HO,46X,21HHYDROGRAPH  FROM  AREA  ,13/) 


n  n 
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is  FORMAT  (1H0,41X, 19HPARTIAL  HYDROGRAPH  ,IU/) 

16  FORMAT  (1HO,39X,29HOUTFLOW  HYDROGRAPH  RESERVOIR  ,U/) 

17  FORMAT  (1HO,A4X,25HOUTFLOW  HYDROGRAPH  REACH  ,U/) 

27  FORMAT(10X,"TIME",6X,"  FLOW" , 1 IX, "TIME" ,6X,"  FLOW" , 1 IX, "TIME" , 
&6X , "FLOW" , 1 IX, "TIME" ,6X , "FLOW" , 1 IX, "TIME" , 6X , "FLOW"/ 1 IX, "HRS" , 
.&7X,"  MS",12X,"HRS",7X,"  MS" , 12X,"HRS" ,7X,"  MS" , 12X,"HRS" , 

&7X,"  MS",12X,"HRS",7X,"  MS") 

19  FORMAT  (5(5X,F10.3,F10.0)) 

21  FORMAT(  'PRINT  HYD' ,T21 ,'ID=' ,11 ,T29,'CODE=' , II) 

26  FORMAT(1HO,9X, "RUNOFF  V0LUME=",F10.0,"  MM  "/10X,"PEAK  DISCHARGE 

&  RATE  =",F10.0,"CMS'7//) 

39  FORMAT  (5(5X,F10.3,F10.2)) 

END 


SUBROUTINE  PUHYD 

THIS  SUBROUTINE  PUNCHES  HYDROGRAPHS  IN  FORM  TO  BE  USED  BY 
SUBROUTINE  RECHD 

COMMON/ BLOCKl /  OCFS ( 300 , 6 ) , DATA( 3 10 ) , CFS ( 300 ) , CTBLE ( 50 , 1 1 ) , 
&RAIN(300),R0IN(6), 

&A(20,6),Q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 
&ZALFA(20) ,IEND(6) ,DA(6) , DIST(6) ,SEGN(6 ) ,DT(6 ) , PEAK(6 ) , ISG(6 ) , 

&NPU , NHD , NER , MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 
DIMENSION  DUMMY (300) 

ID=DATA( 1 ) 

M=IEND(ID) 

IF(ICODE.Eq.O)GO  TO  3 
DA1=DA(ID)*2.590 
PEAK1=PEAK(ID)*0, 02832 
ROINi=ROIN(ID)*25.4 
DO  4  1=1, M 

DUMMY( I )=OCFS( I , ID)*0. 02832 

4  CONTINUE 

WRITE(7,5)ID,NHD,DT(ID),DA1 ,PEAK1,R0IN1 ,IEND(ID) ,ICODE 
WRITE(7,6)(DUMMY(I),I=1 ,M) 

RETURN 

3  WRITE(7, 1)1D,NHD,DT(ID) ,DA(ID) ,PEAK(ID) ,ROIN(ID) ,IEND(ID) ,IC0DE 

WRITE  (7,2)  (0CFS(1,ID),I=1,M) 

RETURN 

C 

I  FORMAT(  'RECALL  HYD' ,T21 ,'ID=' ,11 ,T29,'HYD  NO=' , 13 ,T42 , 'DT=' , F9. 
&6,'  HRS',T61,'DA=',F8.3,'  SQ  MI'/T21 , 'PEAK=' , F7. 0 , 'CFS' ,T40 , 'R0=' , 
5iF6.3,"  INCHES  ",T59,"NO  PTS=" , I3/21X,"CODE=" , II /T2 1 , 

&"FLOW  RATES") 

5  FORMAT(  'RECALL  HYD' ,T2 1 , 'ID=' , II ,T29, 'HYD  NO=' , 13 ,T42 , 'DT=' ,F9. 
&6,'  HRS',T61,'DA=',F8.3,'  SQ  KM'/T2 I , 'PEAK=' , F7 . 2 , 'CMS' ,T40, 'RO=' , 
&F6.0,"  MM  ",T59,"NO  PTS=" , I3/21X, "CODE=" , II /T21 , 


&"FLOW  RATES") 
FORMAT  (T21,7F8.0) 
FORMAT  (T21,7F8.2) 
END 


SUBROUTINE  HPLOT 

THIS  SUBROUTINE  PLOTS  EITHER  1  OR  2  HYDROGRAPHS  ON  A  SET  OF  AXIS 

COMMON /B LOCK 1/  OCFS(300,6) ,DATA(310) ,CFS(300) ,CTBLE( 50 , 1 1) , 
&RAIN(300),R0IN(6), 

&A(20,6),q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 
&ZALFA(20) ,IEND(6) ,DA(6) ,DIST(6) ,SEGN(6) ,DT(6) ,PEAK(6) ,ISG(6) , 
&NPU , NHD , NER ,MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 
ID1=DATA(1) 

ID2=DATA<2) 

DATA  ZERO,  PLUS,  BLANK,  DASH,  DOT/'O' , 

MAX=121 

J=1 

ARE  THERE  1  OR  2  HYDROGRAPHS 
IF  (ID2)  1,1,2 

DETERMINE  HIGHEST  PEAK  IF  2  HYDROGRAPHS 
QMAX=PEAK(ID1) 

GO  TO  14 

IF  (PEAK(ID1)-PEAK(ID2))  3,3,4 
QMAX-PEAK(ID2) 

GO  TO  5 

qMAX=PEAK(IDl) 

IF  2  HYDROGRAPHS  DETERMINE  LARGEST  DT  AND  INTERPOLATE  OTHER 

HYDROGRAPH  IF  NECESSARY 

IF  (DT(I01)-DT(ID2))  6,13,7 

L=ID1 

K=ID2 

GO  TO  8 

L=-ID2 

K=ID1 

M=IEND(L) 

TID=DT(K) 

TIDH=0. 

DO  II  1=2, M 
TIDH=TIDH+DT(L) 

IF  (TID-TIDH)  10,9,11 
J»J+1 

CFS(J)=OCFS(I,L) 

TID*TID+DT(K) 

GO  TO  11 


CFS(J)=0CFS(I-1 ,L)+((TID-TIDH+DT(L))/DT(L))*(0CFS(I,L)-0CFS(I-1 ,L) 

) 

TID=TID+DT(K) 

CONTINUE 

IEND(L)=J 

DT(L)=DT(K) 

DO  12  1=2, J 
OCFS(I,L)=CFS(I) 

IF  (IEND(ID1)-IEND(ID2))  14,14,15 
M=IEND(IDl) 

GO  TO  16 
M=IEND(ID2) 

XM  =  M 

DETERMINE  TIME  SCALE 
XSCL  =  XM  /  120. 

YSCL=QMAX/50. 

PLOT  HYDROGRAPHS 
DO  20  1=1, MAX 
CFS(I)=DASH 
IF(ICODE.EQ.O)GO  TO  49 
WRITE (6, 50) 

F0RMAT(T2,"FL0W  RATE  (CMS)") 

QMAX1=QMAX*0. 02832 

WRITE(6,41  )QM.\X1  ,D0T,(CFS(I),I=1 , MAX), DOT 
GO  TO  51 
WRITE (6, 48) 

F0RMAT(T2,'FL0W  RATE  (CFS)') 
WRITE(6,41)QMAX,D0T,(CFS(I),I=1,MAX),D0T 
Ql-QMAX 
Jl  =  10 

DO  37  J=l,50 
IF  (J-Jl)  23,21,23 
DO  22  1=1, MAX 
CFS(I)=DASH 
GO  TO  25 


DO  24  1=1, MAX 
CFS(I)=BLANK 
Q2=Q1-YSCL 
DO  28  1=2, M 

IF  (OCFS(I,IDl)-qi)  26,27,28 
IF  (OCFS(I,IDl)-q2)  28,28,27 
XI  =  I 

K  -  XI  /  XSCL  +  1. 

CFS(K)=ZLRO 
CONTINUE 

WRITE  (6,44)  DOT,(CFS(I),I=l,MAX),DOT 
IF  (ID2)  34,34,29 
DO  18  I  *  1 ,  MAX 


CFS(I)  =  BLANK 
DO  33  I»1,M 

IF  (OCFS(I,ID2)-Ql)  30,31,33 
IF  (OCFS(I,ID2)-Q2)  33,33,31 
XI  =  I 

K  =  XI  /  XSCL  +  1. 

CFS(K)=PLUS 

CONTINUE 

WRITE  (6,42)  (CFS(I),I=1,MAX) 

IF  (J-Jl)  36,35,36 
J1=J1+10 

IF(ICODE.EQ.O)GO  TO  52 

QD=Q2*0. 02832 

WRITE(6,43)QD 

GO  TO  36 

WRITE(6,43)Q2 

Q1=Q2 

CONTINUE 

CFS(1)=TIME 

DTT=DT(ID1)*(XM  -  1.)  /  12. 

PUT  TIME  ARRAY  IN  CFS  AND  WRITE  TIME  SCALE 

DO  38  1=2,13 

CFS(I)=CFS(I-1)+DTT 

WRITE  (6,45)  (CFS(I),I=1,13) 

WRITE  (6,46) 

IF  (NPU)  40,40,39 
WRITE  (7,47)  IDl,ID2 
RETURN 

FORMAT(1X,F7,0, 123A1) 

F0RMAT(IH+,8X,121A1) 

FORMAT  (1H+,F7.0) 

F0RMAT(8X, 123A1) 

FORMAT(T3, 13F10.2) 

F0RMAT(49X,'TIME  HOURS'///) 

FORMAT/  'PLOT  HYD',T21,'ID  1=' , 1 1 , T29 , ' ID  11  =  ', II) 
END 


SUBROUTINE  ADHYD 

THIS  SUBROUTINE  ADDS  TWO  HYDROGRAPHS. 

COMMON/ BLOCKl/  OCFS( 300 , 6) ,DATA( 310 ) ,CFS( 300 ) ,CTBLE( 50 , 1 1 ) , 
&RAIN(300) ,ROIN(6) , 

&A(20,6),Q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 
&ZALFA(20) ,IEND(6),DA(6),DIST(6),SEGN(6),DT(6) ,PEAK(6),ISG(6), 
&NPU , NHD , NER ,MAXNO , NCOMM, ICC , NCODE .TIME , KCODE , ICODE 


ID=DATA(1) 

NHD=DATA(2) 

ID1=DATA(3) 

ID2=DATA(4) 

PEAK(ID)  =  1. 

MAKE  TIME  INCREMENTS  EQUAL  IF  NOT  EQUAL.  USE  SMALLER  INCREMENT 

IF  (DT(ID1)-DT(ID2))  1,3,2 

DT(ID)=DT(ID1) 

L=ID1 

K=ID2 

GO  TO  6 

DT(ID)=DT(ID2) 

L=ID2 

K=ID1 

GO  TO  6 

DT(ID)=DT(ID1) 

IF  (IEND(ID1)-IEND(ID2))  4,4,5 
M3=IEND(ID1) 

K1=ID2 

IEND(ID)=IEND(ID2) 

GO  TO  18 
M3=IEND(ID2) 

Kl=IDl 

IEND(ID)=IEND(ID1) 

GO  TO  18 

DETERMINE  DUR.4TI0NS  OF  FLOW 
XIEND1=IEND(ID1 )-l 
XIEND2=IEND(ID2)-1 
DUR1=XIEND1*DT(ID1) 

DUR2=XIEND2*DT(ID2) 

IF  (DUR1-DUR2)  7,8,8 
IEND( ID)=DUR2/DT( ID)+1 . 

M3=DUR1/DT(ID)+1. 

K1=ID2 
GO  TO  9 

IEND(ID)=DUR1/DT(ID)+1. 

M3=DUR2/DT(ID)+1. 

Kl=IDi 

IF  (IEND(ID)-300)  11,11,10 

IEND(ID)=300 

M2=IEND(K) 

J=1 

INTERPOLATE  ONE  HYDROGRAPH  IF  NECESSARY 
TIDH=0. 

TID=DT(ID) 

DO  15  1=2, M2 
TIDH=TIDH+DT(K) 

IF  (TIDH-TID)  15,13,14 
J=J+1 

DATA  (J)=OCFS(I,K) 


TID=TID+DT(ID) 

IF  (J-300)  15,16,16 
J=J+1 

DATA  (J)=0CFS(I-1  ,K)+((TID-TIDH+DT(K.))/DT(K))*(0CFS(I,iO-0CFS(I-l 

&K)) 

TID=TID+DT(ID) 

IF  (J-300)  12,16,16 

CONTINUE 

IEND(K)=J 

DO  17  1=2, J 

OCFS(I,K)=DATA(I) 

M=IEND(ID) 

DA(ID)=DA(ID1 )+DA( ID2) 

RO  =  0. 

ADD  HYDROGRAPHS 
DO  20  1=1, M3 

OCFS ( I , ID )=OCFS ( I , ID  1 )+OCFS ( I , ID2 ) 

IF  (OCFS(I,ID)  -  PEAK(ID))  20,20,19 
PEAK(ID)  =  OCFS(I,ID) 

RO  =  RO  +  OCFS(I,ID) 

IF  (PEAK(ID)  -  PEAK(Kl))  21,22,22 
PEAK(ID)  =  PEAK(Kl) 

IF  (M-M3)  25,25,23 
M3  =  M3  +  1 
DO  24  I  =  M3,M 
0CFS(I,ID)  =  OCFS  (I,K1) 

RO  =  RO  +  OCFSd.ID) 

ROIN(ID)  =  (RO  *  DT(ID))  /  (DA(ID)  *  645.333) 

IF  (NPU)  27,27,26 

WRITE  (7,28)  ID,NHD,ID1,ID2 

RETURN 

F0RMAT(  'ADD  HYD' ,T21 , 'ID=' , I 1 ,T29 , '  HYD  r:0=' , 13 ,T45 , ' ID  I=',I1 
&T60,'ID  11=', II) 

END 


SUBROUTINE  SRC 

THIS  SUBROUTINE  STORES  AN  ELEVATION  -  END  AREA  -  FLOW  TABLE. 

COMMON/BLOCKl/  OCFS(300,6) ,DATA(310) ,CFS(300) ,CTBLE( 50 , 1 1 ) , 

5. RAIN(300)  ,ROIN(6)  , 

&A(20,6),Q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 

6. ZALFA(20),IEND(6),DA(6),DIST(6),SEGN(6),DT(6),PEAK(6),ISG(6), 
&NPU , NHD , NER ,MAXNO , NCOMM , ICC , NCODE .TIME , KCODE , ICODE 

ID=0ATA(  1) 

VS=DATA(2) 

VALLEY  SECTION  NUMBER 

REMAINING  DATA  ARE  ELEVATION,  AREA,  AND  FLOW  FOR  EACH  POINT  OF 


THE  RATING  CURVE 
IF(KCODE.EQ.O)GO  TO  2 
J=3 

DO  3  1=1,20 

DATA(J)=DATA(J)/0.3048 

DATA(J+l)=DATA(J+l)/0.093 

DATA(J+2)=DATA(J+2)/0. 02832 

J=J+3 

CONTINUE 

EMIN=DATA(3) 

J=3 

DO  1  1=1,20 

DEEP(I,ID)=DATA(J)-EMIN 

A(I,ID)=DATA(J+1) 

Q(I,ID)=DATA(J+2) 

J=J+3 

CONTINUE 

RETURN 

END 


SUBROUTINE  CMPRC 

THIS  SUBROUTINE  COMPUTES  THE  DISCHARGE  END-AREA  ELEVATION 
RELATIONSHIP  FOR  A  VALLEY  SECTION. 

COMMON/BLOCKl/  OCFS(300,6) ,DATA(310) ,CFS(300) ,CTBLE(50, 1 1 ) , 
&RAIN(300),ROIN(6), 

&A(20,6),q(20,6),DEEP(20,6) ,ITBLE(50,2),DP(20),SCFS(20),C(20) 
&ZALFA(20),IEND(6),DA(6) ,DIST(6),SEGN(6),DT(6),PEAK(6),ISG(6) 
&NPU , NHD , NER ,MAXNO , NCOMM, ICC , NCOOE ,TIME , KCODE , ICODE 
ID=DATA(  1  ) 

STORAGE  LOCATION  NUMBER.  (1-6) 

VS=DATA(2) 

VALLEY  SECTIOI.  IDENTIFICATION  NUMBER. 

NSEG=DATA(3) 

NUMBER  OF  SEGMENTS  IN  THE  VALLEY  SECTION. 

IF ( KCODE. EQ.O)GO  TO  26 
DATA(4)=DATA(4)/0.3048 
DATA(5)=DATA(5)/0.3048 
ELO=DATA(4) 

EMAX=DATA(5) 

MAXIMUM  ELEVATION  FOR  COMPUTATIONS. 

SL0PE1=DATA(6) 

CHANNEL  SLOPE. 

SLOPE2=DATA(7) 

FLOOD  PLAIN  SLOPE. 

DIF=(EMAX-ELO)/19. 
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ec  1)=EL0 
DO  1  1=2,20 

1  C(I)=C(I-1)+DIF 

C  SET  AREA  AND  DISCHARGE  ARRAYS  =  0. 

DO  2  1=1,20 

A(r,rD)=o. 

2  Q(I,ID)=0. 

J=8 

C  READ  N  VALUES  AND  SEGMENT  BORDER  POINTS. 

DO  3  I=1,NSEG 
SEGN(I)=DATA(J) 

IF(KCODE.NE.0)DATA(J+l)=DATA(J+l)/0.3048 

DIST(I)=DATA(J+1) 

3  J=J+2 

C  REMAINING  DATA  ITEMS  ARE  DISTANCES  AND  ELEVATIONS 

IF(KCODE.EQ.O)GO  TO  27 
DO  28  I=J,310 
DATA(I)=DATA(I)/0.3048 
28  CONTINUE 

27  JJJ=J 

DO  6  I=1,NSEG 

4  J=J+2 

IF  (DATACJ)  -  DIST(I))  4,5,5 

5  ISG(I)  =  J  +  1 

6  CONTINUE 

C  COMPUTE  DISCHARGES  AND  END  AREAS  FOR  EACH  SEGMENT 

DO  22  K«1,NSEG 
J=JJJ 
JJJ1=JJJ+1 
IF  (SEGN(K))  7,7,8 

7  SL0PE=SL0PE1 
SEGN(K)=-SEGN(K) 

GO  TO  9 

8  SLOPE=SLOPE2 

9  SLPN=1.486*SL0PE**.5 

C  COMPUTE  AREA  AND  DISCHARGE  FOR  SEGMENT. 

DO  21  1=2,20 
AA=0. 

P=0. 

J=JJJ-1 

DEP2=0. 

10  J=J+2 

IF  (J-ISG(K))  12,12,11 

11  IF  (AA-.OOl)  21,21,20 

12  IF  (DATA(J)-C(I))  13,10,10 

13  DEP1=C(I)-DATA(J) 

IF  (J-JJJl)  16,16,14 

14  XL-DATA(J-l)-DATA(J-3) 

DE  P3  =ABS ( DATA ( J -2 ) -DATA ( J ) ) 

XL=XL*DEP1/DEP3 


AA=AA+XL* ( DE P 1 +DE P 2 ) / 2 . 

P=P+SQRT ( ( DEP 1 -DEP2 ) **2+XL**2 ) 

DEP2=DEP1 

J=J+2 

IF  (J-ISG(K))  17,17,20 
IF  (DATA(J)-C(I))  18,18,19 
DEP1-C(I)-DATA(J) 

XL=DATA(J-1 )-DATA(J-3) 

GO  TO  15 
DEP 1=0. 

XL=DATA(J-l)-DATA(J-3) 

DEP3=ABS ( DATA( J-2 )-DATA( J ) ) 

XL=XL*DEP2/DEP3 
AA=AA+XL* ( DE  P 1 +DE  P2 ) / 2 . 

P=P+SQRT((DEP1-DEP2)**2+XL**2) 

DEP2=0. 

GO  TO  10 
R=AA/P 

SGN=SEGN(K)  -  .0025*R 

ADD  DISCHARGES  AND  AREAS  FOR  ALL  SEGMENTS  TO  OBTAIN  TOTALS  FOR 
VALLEY  SECTION. 

Q( I , ID )=Q( I , ID)+AA*R** . 66667*SLPN/SGN 

A(I,ID)=A(I,ID)+AA 

CONTINUE 

JJJ=J-3 

CONTINUE 

IF(ICODE.EQ.O)GO  TO  29 

WRITE(6,31)VS 

DO  30  1=1,20 

C1=C(I)*0.3048 

A1=A(I,ID)*0.093 

qi=Q(I,ID)*0. 02832 

DEEP(I,ID)=C(I)-ELO 

WRITE(6,32)Cl,Al,qi 

CONTINUE 

RETURN 

WRITE(6,24)VS 

DO  23  1=1,20 

DEEP(I,ID)=C(I)-ELO 

WRITE  (6,25)  C(I),A(I,ID),q(I,ID) 

CONTINUE 

RETURN 

FORMAT(1HO,T42, 'RATING  CURVE  VALLEY  SECTION' , F5 . I /T46 , 'WATER' ,T56 , 
&'FLOW', T66, 'FLOW' /T45, 'SURF ACE', T56,' AREA' ,T66 , 'RATE' /T46 , 'ELEV' , 
&T56,'Sq  FT',T66,'CFS') 

FORMAT (1 HO, T42,' RATING  CURVE  VALLEY  SECTION' , F5 . I /T46 , 'WATER' , T56 , 
&'FL0W',T66,'FL0W'/T45,'SURFACE',T56,'AR£A',T«36, 'RATE' 'T^h. 'ELEV' , 
&T56,'Sq  M' ,T66,'CMS') 

FORMAT  (40X,F10.2,2F10. 1 ) 


FORMAT  (40X,3F10.2) 
END 


SUBROUTINE  STT 

THIS  SUBROUTINE  STORES  A  DEPTH  -  FLOW  -  TRAVEL  TIME  TABLE. 

COMMON/BLOCKl /  OCFS ( 300 , 6 ) , DATA( 3 1 0 ) , CFS ( 300 ) , CTBLE ( 50 , 1 1 ) , 
&RAIN(300),ROIN(6), 

&A(20,6),Q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20) 
&ZALFA(20) ,IEND(6) ,DA(6) ,DIST(6) ,SEGN(6) ,DT(6) ,PEAK(6) ,ISG(6) 
&NPU , NHD , NER , MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 
ID*DATA( I ) 

REACH=DATA(2) 

MET1=DATA(5) 

IF(METl.EQ.O)GO  TO  2 

DATA(3)=DATA(3)/0.3048 

J=6 

DO  3  1=1,19 

DATA(J)-DATA(J)/0.3048 
DATA(J+l)=DATA(J+l)/0. 02832 
J=J+3 

XL=DATA(3) 

SL0PE=DATA(4) 

DIST(ID)=SLOPE*XL 

J=6 

DO  1  1=1,19 
DP(I)-DATA(J) 

SCFSU)=DATA(J+1) 

C(I)=DATA(J+2) 

J=J+3 

RETURN 

END 


SUBROUTINE  CMPTT 

THIS  SUBROUTINE  COMPUTES  THE  TRAVEL  TIME  AT  GIVEN 
DISCHARGE  RATES 

COMMON/BLOCKl /  OCFS ( 300 , 6 ) , DATA( 3 10 ) , CFS ( 300 ) , CTBLE ( 50 , 1 1 ) , 
&RAIN(300) ,ROIN(6) , 

&A(20,6),Q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20) 
&ZALFA(20),IEND(6),DA(6),DIST(6),SEGN(6),DT(6),PEAK(6),ISG(6) 
&NPU , NHD , NER , MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 
ID=DATA( 1 ) 

REACH=DATA(2) 

NOVS=DATA(3) 


IF ( KCODE . NE . 0 ) DATA( 4 ) -DATA( 4 ) / 0 . 3048 
XL-DATA(4) 

SL0PE=DATA(5) 

DIST(ID)=SLOPE*XL 
XLD36  -  XL  /  3600. 

ZERO  ARRAYS 
DO  1  J=l,20 
DATA  (J)=0. 

CFS(J)-0. 

ID1=1 

FIND  RATING  CURVE  WITH  SMALLEST  MAXIMUM  FLOW  RATE 
QMIN=Q(20,ID1) 

MIN=ID1 
GO  TO  4 
IDl-IDl+1 

IF  (QMIN-Q(20,ID1))  4,4,2 
IF  (IDl-NOVS)  3,5,5 
I-l 

SET  SCFS  ARRAY  EQUAL  TO  Q  ARRAY  OF  LOWEST  RATING  CURVE 

DO  6  J=2,20 

SCFS(I)=Q(J,MIN) 

1=1+1 

COMPUT  END  AREA  AND  DEPTH 
DO  9  IDl=l,NOVS 
DO  9  J=l,19 
DO  7  1=2,20 

IF  (Q(I,IDl)-SCFS(J))  7,17,8 
CONTINUE 

DATA  (J)=A(I,ID1)+DATA(J) 

CFS(J)=DEEP(I,ID1)+CFS(J) 

GO  TO  9 

XY=(SCFS(J)-q(I-l,IDl))/(Q(I,IDl)-Q(I-l,IDl)) 

DATA  ( J )=A( I-l , ID  I )+XY*( A( 1,101 )-A( I-l , IDl) )+DATA( J) 

CFS( J)=DEEP( I-l , IDl )+XY*(DEEP( I , IDl )-DEEP(I-l , IDl ) )+CFS( J) 

CONTINUE 

XNOVS=NOVS 

IF(ICODE.EQ.O)GO  TO  19 
WRITE (6, 20) REACH 
GO  TO  21 

WRITE(6,13)REACH 
DO  10  1=1,19 

AVAREA  =  DATA  (I)  /  XNOVS 
DP  (I)  =  CFS(I)  /  XNOVS 
S  =  AVAREA  *  XLD36 
C(I)  =  S/SCFS(I) 

IF(ICODE.EQ.O)GO  TO  24 
DP1=DP(I)*0.3048 
SCFS1=SCFS(I)*0. 02832 
WRITE(6,14)DP1,SCFS1,C(I) 


GO  TO  10 

WRITE(6,14)DP(I),SCFS(I),C(I) 

CONTINUE 
PUNCH  CODE 
IF(NPU)12, 12,25 
IF(ICODE.EQ.O)GO  TO  11 
XL1-XL*0.3048 

WRITE (7 , 22 ) ID , REACH , XLl .SLOPE , ICODE 
DO  23  1-1,19 
DP1=DP(I)*0.3048 
SCFS1-SCFS(I)*0. 02832 
WRITE(7,26)DP1,SCFS1,C(I) 

CONTINUE 

RETURN 

WRITE ( 7 , 1 5 ) ID , REACH , XL , SLOPE , ICODE 
WRITE  (7,16)  (DP(I),SCFS(I),C(I), 1-1,19) 

RETURN 

FORMAT(1HO,T46, 'TRAVEL  TIME  TABLE'/T54, 'REACH' ,F5. 1//T46, 'WATER' ,T 
&56,'FLOW' ,T65,'TRAVEL'/T46, 'DEPTH', T56, 'RATE', T66, 'TIME' /T46, 'FEET 
&',T56,'CFS',T66,'HRS') 

FORMAT  (40X,F10.2,F10.0,F10.2) 

FORMAT('STORE  TRAVEL  TIME' ,T2l , 'ID-' , II , T29 , 'REACH  NO-' ,F5. 1 ,T44 , 
&'LENGTH-' ,F9.0,'  FT'/T21 , 'SLOPE-' , F8. 6 , 'FT/FT' , "CODE-" , II /T2 
&1,'DEPTH(FT)' ,T35,'FLOW(CFS)',T49,'TIME(HRS)') 

F0RMAT(1H0,T46, 'TRAVEL  TIME  TABLE'/T54 , 'REACH' , F5. 1//T46 , 'WATER' ,T 
456 , 'FLOW' ,T65 , 'TRAVEL'/T46, 'DEPTH' ,T56, 'RATE' , T66 , 'TIME'/T46 , 
&"METER" ,T56 , 'CMS' ,T66 , 'HRS' ) 

FORMATC'STORE  TRAVEL  TIME' ,T2l , ' ID-' , I 1 ,T29 , 'REACH  NO-' ,F5. I ,T44 , 
4' LENGTH-', F9.0,'  M'/T21 , 'SLOPE-' , F8. 6 , 'M/M' ."CODE-" , 1 1 /T2 
4l,'DEPTH(M)',T35,'FL0W(CMS)',T49,'TIME(HRS)') 

FORMAT  (T21 ,F7.2,F15.2,F15.3) 

FORMAT  (T21 ,F7,2,2F15.3) 

END 


SUBROUTINE  ROUTE 

THIS  SUBROUTINE  ROUTES  A  HYDROGRAPH  THROUGH  A  REACH  WITH  THE 
NEW  VSC  METHOD  OF  FLOOD  ROUTING.  THIS  METHOD  ACCOUNTS  FOR  THE 
VARIATION  IN  WATER  SURFACE  SLOPE. 

COMMON/ BLOCKl /  OCFS ( 300 , 6 ) , DATA( 3 1 0 ) , CFS ( 300 ) , CTBLE ( 50 , 1 1 ) , 
4RAIN(300) ,ROIN(6) , 

4A(20,6),Q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 
4ZALFA( 20 ) , IEND( 6 ) , DA( 6 ) , DIST( 6 ) , SEGN( 6 ) , DT ( 6 ) , PEAK( 6 ) , ISG( 6 ) , 
4NPU , NHD , NER ,MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 
ID-DATA ( I ) 

NHD=DATA(2) 

IDH-DATA(3) 


DT(ID)»DATA(4) 

DA(ID)=DA(IDH) 

M=IEND(IDH) 

IF  ID  AND  IDH  ARE  EQUAL,  ADD  1  TO  IDH 
IF  (ID-IDH)  3,1,3 
IDH=IDH+1 
DO  2  1=1, M 

OCFS(I,IDH)-OCFS(I, IDH-1 ) 

DT(IDH)-DT(IDH-l ) 

PEAK(IDH)=PEAK( IDH-1  ) 

NERRT-0 
PEAK(ID)  -  1. 

RO  -  0. 

N-19 

OCFS(1,ID)-0. 

S  -  0. 

T1  -  C(l) 

J«1 

CUES  -  1. 

CFS(l)-0. 

IF  ROUTING  INTERVAL  IS  NOT  EQUAL  TO  TIME  INCREMENT  OF  INFLOW 
HYDROCRAPH,  INTERPOLATE 
IF  (DT(ID)-DT(IDH))  8,15,4 
TID-DT( ID) 

TIDH-0. 

DO  7  1-2, M 
TIDH-TIDH-t-DT(IDH') 

IF  (TID-TIOH)  6,5,7 
J-J  +  l 

CFS(J)-OCFS(I,IDH) 

TID-TID+DT( ID) 

GO  TO  7 
J-J+1 

CFS(J)=OCFS(I-l ,IDH)+((TID-TIDH+DT(IDH))/DT(IDH))*(OCFS(I,IDH)-OCF 
4S( I-l , IDH)) 

TID=TID+DT( ID) 

CONTINUE 
GO  TO  13 
TIDH-0. 

TID-DT(ID) 

DO  12  1-2, M 
TIDH-TIDH+DTdDH) 

IF  (TIDH-TID)  12,10,11 
J-J+1 

CFS(J)-OCFS(I,IDH) 

TID-TID+DT(ID) 

IF  (J-300)  12,13,13 
J-J+1 

CFS(J)=OCFS(I-l ,IDH)+((TID-TIDH+DT(IDH))/DT(IDH))*(OCFS(I,IDH)-OCF 
SiS(I-l,IDH)) 
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TID-TID+DT(ID) 

IF  (J-300)  9,13,13 

12  CONTINUE 

13  IEND(IDH)-J 
DT(IDH)-DT(ID) 

M-J 

DO  14  1-2, M 

14  OCFS(I,IDH)-CFS(I) 

C  IF  INFLOW  IS  ZERO,  SO  IS  OUTFLOW 

15  DO  16  L-2,M 

IF  (OCFS(L,IDH))  16,16,49 

16  0CFS(L,ID)-0. 

C  ROUTE 

49  DATA  (L-1)  -  0. 

DO  42  I-L,300 
IF  (I-M)  18,18,17 

17  OCFS(I,IDH)-OCFS(I-l ,IDH)*.9 

18  AVIN-(OCFS(I,IDH)+OCFS(I-l,IDH))/2. 

SIA  -  AVIN  +  S 

J-1 

C  DETERMINE  DEPTH  AND  TRAVEL  TIME  OF  INFLOW 

IF  (OCFS(I,IDH)-SCFS(l))  19,23,20 

19  DI2  -  (0CFS(I,IDH)  /  SCFS(l))  *  DP(1) 

TI2  -  C(l) 

GO  TO  25 

20  DO  21  J=2,N 

IF  (OCFS(I,IDH)-SCFS(J))  24,23,21 

21  CONTINUE 

IF  (NERRT)  22,22,36 

22  WRITE  (6,46) 

NERRT- 1 

GO  TO  36 

23  DI2-DP(J) 

TI2  -  C(J) 

GO  TO  25 

24  RATIO=(OCFS(I,IDH)-SCFS(J-1))/(SCFS(J)-SCFS(J-1 )) 
DI2-DP(J-l)+RATIO*(DP(J)-DP(J-l)) 
TI2=C(J-l)+RATIO*(C(J)-C(J-l)) 

25  DO  35  IT-1,10 
J=1 

C  DETERMINE  DEPTH  AND  TRAVEL  TIME  OF  OUTFLOW 

IF  (GUES-SCFS(l))  26,29,27 

26  D02  =  (GUES  /  SCFS(l))*  DP(1) 

T02  =  C(l) 

GO  TO  31 

27  DO  28  J=2,N 

IF  (GUES-SCFS(J))  30,29,28 

28  CONTINUE 
J-N 

29  D02-DP(J) 


T02=C(J)  i 

GO  TO  31 

RATIO=(GUES-SCFS ( J-1 ) ) / ( SCFS ( J )-SCFS ( J-1 ) ) 

IX)2=DP(  J-1  )+RATIO*(DP(  J)-DP(  J-D) 

T02=C(J-1 )+RATIO*(C(J)-C( J-1 )) 

FIND  WATER  SURFACE  SLOPE 
DDD-DIST(ID)/(DIST(ID)+DI2-D02) 

IF  (DDD-.Ol)  32,32,33 
GUES=0CFS(I-1,IDH) 

GO  TO  35 

T2  =  .5  *  (TI2  +  T02) 

T2=T2*SQRT(DDD) 

T  -  T1  +  T2 

COMPUTE  ROUTING  COEFFICIENT 
COEF  =(2.  *  DT(ID))  /  (T+DT(ID)) 

02  -  COEF  *  SIA 
TRYl  =  CUES 
RATI0=02/(GUES+. lE-20) 

DIFF=ABS(1. -RATIO) 

TEST  FOR  CONVERGENCE 
IF  (DIFF-.OOl)  37,37,34 
GUES-02 
CONTINUE 

0CFS(I,ID)=-DATA(I-1)*SIA 
DATA(I)  =  DATA(I-l) 

WRITE  (6,47)  I,OCFS(I,ID) 

GO  TO  38 

0CFS(I,ID)=DATA(I-1)*SIA 
DATA(I)  »  DATA(I-l) 

GO  TO  38 
0CFS(I,ID)=02 
DATA  (I)  =  COEF 
COMPUTE  NEW  STORAGE 
S  =  SIA  -  OCFS(I,ID) 

T1  =  T2 

RO  =  RO  +  OCFS  (I, ID) 

IF  (OCFS(I,ID)  -  0CFS(I-1,ID))  39,40,40 
IF(OCFS(l,ID)  -1.)  43,43,42 
IF(OCFS(i,ID)  -  PEAK(ID))  42,42,41 
PEAK(ID)=OCFS  (I, ID) 

CONTINUE 


1-300 

IEND(ID)=«I 

ROIN(ID)  -  (RO*DT(ID))/(DA(ID)*645.333) 

PUNCH  CODE 

IF  (NPU)  45,45,44 

WRITE  (7,48)  ID,NHD,IDH,DT(ID) 

RETURN 


o  o 


46  FORMATdHO,  'TRAVEL  TIME  TABLE  EXCEEDED') 

47  FORMAT(T10, 'PROBLEM  FAILED  TO  CONVERGE  AFTERIO  ITERATIONS.  CONVERG 
&ENCE  WAS  FORCED. '/T20, 'OUTFLOW  NUMBER  =  ',1^,'RATE  =',F10.2) 

48  FORMAT(  'ROUTE' ,T21 ID=' , II ,T29 ,'HYD  NO=' , 13 ,T45 ,' INFLOW  ID=',I 
&1,T65,'DT=',F8.6,'HRS') 

END 


SUBROUTINE  RESVO 

THIS  SUBROUTINE  ROUTES  A  HYDROGRAPH  THROUGH  A  RESERVOIR  WITH  THE 
STORAGE-INDICATION  METHOD. 

COMMON/BLOCKl /  OCFS (300,6), DATA( 31 0) , CFS ( 300 ) , CTBLE (50,11), 
&RAIN(300),ROIN(6) , 

&A(20,6),Q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 
&ZALFA(20),IEND(6) ,DA(6),DIST(6),SEGN(6),DT(6),PEAK(6),ISG(6), 
&NPU , NHD , NER , MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 
ID=DATA(1) 

NHD=DATA(2) 

IDH=DATA(3) 

NERES=0 

DT(ID)=DT(IDH) 

RO  =  0. 

DA(ID)=DA(IDH) 

PEAK(ID)  =  1. 

J=l 
1=4 

C  REMAINING  DATA  ARE  FLOW  AND  STORAGE  VALUES 

IF(KCODE.EQ.O)GO  TO  25 
DATA(I)=DATA(I)/0. 02832 
DATA(I+1)=DATA(I+1)/1.2I968 

25  SCFS(J)=DATA(I) 

ST0RE1=DATA(I+1)*12. I 
ST0RE=ST0RE1 

C  COMPUTE  STORAGE  COEFFICIENT  ARRAY  C 

1  C(J)=(SCFS(J)/2.)+(ST0RE/DT(ID)) 

1=1+2 

J=J+1 

IF  (J-20)  2,2,3 

2  IF(KCODE.EQ.O)GO  TO  26 
DATA(I)=DATA(I)/0.02832 
DATA(I+l)=DATA(I+l)/ 1.2 1968 

26  SCFS(J)=DATA(I) 

STOR£=DATA(I+l)*12. 1 

IF  (SCFS(J)-.OOl)  3,3,1 

3  N=J-1 
OCFS(1,ID)=0. 

S=STOREl/DT(ID) 
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C  ROUTE 

DO  15  1=2,150 
IF  (I-IEND(IDH))  5,5,4 

4  OCFS(I,IDH)=0.0 

5  AVIN=(0CFS(I,IDH)+0CFS(l-l,IDH))/2. 

SIA=S+AVIN 

C  DETERMINE  PROPER  C 

DO  6  J=1,N 
IF  (SIA-C(J))  10,9,6 

6  CONTINUE 

IF  (NERES)  7,7,8 

7  WRITE  (6,19) 

NERES=1 

8  RESC=SCFS(N)/C(N) 

C  COMPUT  OUTFLOW 

OCFS(I,ID)=RESC*SIA 
GO  TO  11 

9  OCFS(I,ID)=SCFS(J) 

GO  TO  11 

10  OCFS(I,ID)=SCFS(J-l)+((SIA-C(J-l))/(C(J)-C(J-l)))*(SCFS(J)-SCFS(J- 

&1)) 

C  DETERMINE  NEW  STORAGE 

11  S=SIA-OCFS(I,ID) 

RO  =  RO  +  OCFS(I,ID) 

IF  (0CFS(I,ID)-0CFS(I-1,ID))  12,13,13 

12  IF  (0CFS(I,ID)-1.)  16,16,15 

13  IF(OCFS(I,ID)  -  PEAK(ID))  15,15,14 

14  PEAK(ID)  =  OCFS(I,ID) 

15  CONTINUE 
1=150 

16  IEND(ID)=I 

ROIN(ID)  =  RO  *  DT(ID)/(DA(ID)*645.333) 

C  PUNCH  CODE 

IF  (NPU)  18,18,17 

17  II=2*N+3 
IF(ICODE.EQ.O)GO  TO  22 
WRITE ( 7 , 24 ) ID , NHD , IDH , KCODE 
DO  23  1=5,11,2 
DATA(I)=DATA(I)*0.02832 
DATA(I+1 )=DATA(I+1 )* 1.21968 

23  CONTINUE 

WRITE(7,27)(DATA(I),I=5,II) 

RETURN 

22  WRITE (7, 20) ID, NHD, IDH, ICODE 

WRITE  (7,21)  (DATA(I),I=5,II) 

18  RETURN 
C 

19  FORMAT  (1HO,33HSTORAGE-DISCHARGE  TABLE  EXCEEDED.) 

20  FORMAT(  'ROUTE  RESERVOIR' ,T2 1 , 'ID=' , 1 1 ,T29 , 'HYD  NO=' , 13 ,T42 , ' INF 

&LOW  ID=',I1,T60,"CODE=",I1  /T21 , 'OUTFLOW(CFS)' ,T37 , 'STOR 


I 


&AGE(AC  FT)') 

24  FORMAT(  "ROUTE  RESERVOIR' ,T21 , 'ID=' , II ,T29 , 'HYD  NO=' , 13 ,T42 , 'INF 
&LOW  ID=' ,11 ,T60,"CODE='',I1  /T21 ,'OUTFLOW(CMS)' ,T37 ,'STOR 

&AGE(lOOOCU  M)') 

21  FORMAT  (T21.F10.1,F13.1) 

27  FORMAT  (T21 ,F10.2,F13.2) 

END 


SUBROUTINE  ERROR 

C  This  subroutine  determines  the  error  standard  deviation  and  the  peak  flow 
C  error  for  2  hydrographs  (original  program  retained). 

C  Assumes  that  measured  is  IDl 

C  In  addition,  10  other  measures  of  goodness  of  fit  are  calculated. 

C  All  indicies  are  printed  out  in  metric  units. 

COMMON/ BLOCKl/  OCFS ( 300 ,6) ,DATA(310) ,CFS(300) ,CTBLE( 50 , 1 1 ) , 
&RAIN(300),ROIN(6), 

&A(20,6),Q(20,6) ,DEEP(20,6),ITBlE(50,2) ,DP(20) ,SCFS(20) ,C(20) , 
&ZALFA(20) ,IEND(6) ,DA(6) ,DIST(6) ,SEGN(6),DT(6) ,PEAK(6) ,ISG(6) , 

&NPU , NHD , NER , MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 

ID1=DATA(1  ) 

ID2=DATA(2) 

SSE=0. 

WRITE(6,21) 

21  F0RMAT( IH0,T33, 'TIME', T55, 'FLOW  l',T76, 

&  'FLOW  2' ,T95,'ERR0R'/T34, 

&  'HRS' ,T57 ,'CMS' ,T78,'CMS' ,T97 ,'CMS') 

22  J=1 

C  If  the  time  increments  are  not  equal,  interpolate. 

IF  (DT(ID1)-DT(ID2))  1,8,2 

1  L=IDI 
iC=ID2 
GO  TO  3 

2  L=ID2 
K=ID1 

3  M=IEND(L) 

TID=DT(K) 

TIDH=0. 

DO  6  1=2, M 
TIDH=TIDH+DT(L) 

IF  (TID-TIDH)  5,4,6 

4  J=J+1 
CFS(J)=OCFS(I,L) 

TID=TID+DT(K) 

GO  TO  6 


5 


J=J+1 

CFS(J)=0CFS(I-1 ,L)+((TID-TIDH+DT(L))/DT(L))*(0CFS(I,L)-0CFS(I-1 ,L) 

&) 

TID=TID+DT(K) 

6  CONTINUE 
IEND(L)=J 
DT(L)=DT(K) 

DO  7  1=2, J 

7  OCFS(I,L)=CFS(I) 

8  IF  (IEND(ID1)-IEND(ID2))  9,9,10 

9  M=IEND(ID1) 

GO  TO  11 

10  M=IEND(ID2) 

11  T2=TIME 

IF  (ICCODE.EQ.O)THEN 
DO  997  1=1, M 

0CFS(I,ID1)=0CFS(I,ID1)*.02832 
997  OCFS(I,ID2)=OCFS(I,ID2)*. 02832 

ENDIF 


C  Determine  error  -  original  method 
DO  12  1=1, M 

ERR=0CFS(I,ID1)-0CFS(I,ID2) 

WRITE(6,16)T2,OCFS(I,IDl),OCFS(I,ID2) ,ERR 
16  FORMAT  (6X,F12.3,3F12.0) 

25  T2=T2+DT(ID1) 


C  Sum  of  squares  of  error 

12  SSE=SSE+ERR*ERR 
XM=M 


C  Error  variance 
EVAR=SSE/XM 

C  Error  standard  deviation 
ESDEV=SqRT(EVAR) 

C  Percent  error  for  peak  discharge 

ERPK=ABS(PEAK(IDl)-PEAK( ID2)) 
PCTER=(ERPK/PEAK(ID1))*100. 

C  Other  goodness  of  fit  calculations... 


SUMO 1=0. 


SUM0=0. 

SUM1=0. 

SUM2=0. 

SUM3=0. 

SUM4=0. 

SUM5=0. 

SUM6=0. 

SirM7=0. 

SUM8=0. 

SUM9=0. 

SUM10=0. 

SUM11=0. 

SUM12=0. 

DO  77  1=1, M 

ERR=OCFS ( I , ID  1 )-OCFS ( I ,  ID2  ) 

rF(OCFS(I,rDl).EQ.O.O.AND.OCFS(I,ID2).NE.O.O)THEN 
LOGERR=ALOG (OCFS(I,ID2)) 

ELSE  IF(OCFS(I,ID1).NE.O.O.AND.OCFS(I,ID2).EQ.O.O)THEN 
LOGERR=ALOG (OCFS(I,IDl)) 

ELSE  IF(OCFS(I,ID1).EQ.O.O.AND.OCFS(I,ID2).EQ.O.O)THEN 
LOGERR=0. 

ELSE 

LOGERR=ALOG(OCFS(I , IDl ) )-ALOG(OCFS ( I , ID2) ) 

ENDIF 

SUMO=OCFS ( I , ID  1 )+SUMO 

SUM01=OCFS(I,ID2)+SUM01 

SUML=ERR+SUM1 

SUM2=ERR**2+SUM2 

SUM3=LOGERR**2+SUM3 

IF(OCFS(I,IDl),EQ,0.)OCFS(I,IDl)=l. 

SUM4=( ( ERR/OCFS ( I , ID  1 ) )**2 )+SUM4 
CONTINUE 

DO  13  1=2, M 

DIFFL=OCFS(I,IDl)-OCFS(I-I,IDl) 

DIFF2=OCFS(I , ID2)-OCFS( I-l , ID2) 

SUM5=((DIFFl-DIFF2)**2)+SUM5 

SUM7=DIFF1+SUM7 

IF(DIFF1.EQ.0.)DIFF1=1. 

SUM6=(((DIFF1-DIFF2)/DIFF1)**2)+SUM6 

CONTINUE 


SIMMEAN=SUM01/M 

OBSMEAN=SUMO/M 

DIFFMI=SUM7/M 


DO  14  1=2, M 

SUM8=(((0CFS(I,ID1)-0CFS(I-1 ,ID1))-DIFFM1 )**2)+SUM8 


SUM9=((((0CFS(I,ID1)-0CFS(I-1,ID1))/DIFFM1)-1)**2)+SUM9 

CONTINUE 

DO  73  1=1, M 

SUM10=((0CFS(I,ID1)-0BSMEAN)**2)+SUM10 
SUMl 1  =  ( ( (OCFS ( I , ID  1) /OBSMEAN)-! )**2 )+SUMl 1 
SUM12=((OCFS(I,ID2)-SIMMEAN)**2)+SUM12 
CONTINUE 

SDM=SQRT(SUM10/(M-1)) 

SDS=SQRT(SUM12/(M-1)) 

DO  115  1=1, M 

SUM13=( (OCFSd , IDl )-OBSMEAN)/SDM)*( (OCFS( I , ID2)- 
SIMMEAN)/SDS)+SUM13 


0F1=SUM1 

OF2=SUM2 

OF3=SUM3 

OF4=SUM4 

OF5=SUM5 

OF6=SUM6 

OF7=SUM2/SUM10 

OF8=SUM4/SUMll 

OF9=SUM5/SUM8 

OF10=SUM6/SUM9 

AM=M 

OFll-(l./AM)*SUM13 
WRITE (6, 95) 

FORMAT( IHO, lOX,' - 

WRITE (6, 50) 

FORMAT(15X,'  MEASURES  OF  FIT  '//) 

WRITE(6,91) 

FORMATC  lOX,' - ') 

WRITE(6,51 )0F1 


FORMATC 1 OX, 'SUM  OF  ERRORS 
WRITE(6,52)OF2 

',F20.5) 

FORMATC 10X,'OLSQ 

WRITEC6,53)OF3 

',F20.5) 

FORMATC 1 OX, 'LOG  LSQ 
WRITEC6,54)OF4 

',F20.5) 

FORMATC lOX, 'RELATIVE  ERROR 
WRITEC6,55)OF5 

',F20.5) 

FORMATC 1 OX,' ABS  ERROR  -  DIFF 
WRITEC6,56)OF6 

',F20.5) 

FORMATC 1 OX,' RE L  ERROR  -  DIFF 
WRITEC6,57)OF7 

',F20.5) 

FORMATC 1 OX, 'ABS  ERROR/VAR 

',F20.5) 
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WRITE  (6,58  )0F8 

58  FORMAT ( 1 OX,' RE L  ERROR/VAR  ',F20.5) 

WRITE(6,59)OF9 

59  FORMAT(10X,'ABS  ERROR(dif f )/VAR  ',F20.5) 

WRITE(6,60)OF10 

60  FORMAT (1 OX, 'REL  ERROR(dif f )/VAR  ',F20.5) 

WRITE(6,61)0F11 

61  ■ FORMAT (1 OX, 'PEARSONS  r  ',F20.5) 

WRITE(6,92)ESDEV 

92  FORMAT(10X,'ERR  STANDARD  DEV  ',F20.5) 

WRITE(6,93)PCTER 

93  FORMAT(10X,'PEAK  Q  ERROR  ',F20.5) 

WRITE(6,96) 

96  FORMATdOX,' - ') 


WRITE  (6,98) 

98  FORMAT  (//10X,'NOTE:  All  indicies  are  in  metric  units') 

IF  (KCODE.EQ.O)THEN 
DO  9969  1=1, M 

OCFS(I,ID1)=OCFS(I,ID1)/.02832 
9969  OCFS(I,ID2)=OCFS(I,ID2)/.02832 

ENDIF 

RETURN 

C 

END 


SUBROUTINE  SEDT 

C  THIS  SUBROUTINE  COMPUTES  THE  SEDIMENT  YIELD  FOR  A  FLOOD 

COMMON/ B LOCKl /  OCFS ( 300 , 6 ) , DATA( 3 10 ) , CFS ( 300 ) , CTBLE ( 50 , 11 ) , 
&RAIN(300),ROIN(6), 

&A(20,6),Q(20,6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 
&ZALFA(20) ,IEND(6) ,DA(6),DIST(6) ,SEGN(6) ,DT(6) ,PEAK(6) ,ISG(6) , 
&NPU , NHD , NER , MAXNO , NCOMM, ICC , NCODE , TIME , KCODE , ICODE 
ID=DATA(1) 

S0IL=DATA(2) 

CR0P=DATA(3) 

CP=DATA(4) 

SL=DATA(5) 

C  COMPUTE  SEDIMENT  YIELD 

X=ROIN(ID)*DA(ID)*53.333*PEAK(ID) 

SED=95 . *X** . 56*SO  rL*CROP*CP*SL 
IF(ICODE.EQ.O)GO  TO  5 
SED1=SED*0.9072 
WRITE(6,6)SED1 
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GO  TO  7 

5  WRITE  (6,3)  SED 

C  PUNCH  CODE 

7  IF(NPU)2,2,1 

1  WRITE  (7,4)  ID, SOIL, CROP, CP, SL 

2  RETURN 

3  FORMAT  (lOX,  'SEDIMENT  YIELD  -  ',  FlO.l,  '  TONS*' 

4  FORMAT(  'SEDIMENT  YIELD' ,T2 1  ,  '  ID-' ,  1 1  ,T29  ,  '  SO  1 1.- '  .  F5  .  i ,  T-*:  ,  '  R' )P 
&»',F5.3,T57,'CP-',F5.3,T70,'LS-',F5. 3) 

6  FORMATdOX, "SEDIMENT  YIELD-", FIO.  1  ."METRIC  TON") 

END 


BLOCK  DATA 

C  BLOCK  DATA  SUBPROGRAM  UZED  TO  INITIALIZE  ZALFA.CTBLE , ITBLE 
C  AND  NCOMM. 

COMMON/ BLOCKl/  OCFS(300 ,6) ,DATA( 310) ,CFS(300) ,CTBLE( 50 , 1 1 ) , 

&RAIN(300),ROIN(6),  I 

&A(20,6),Q(20.6),DEEP(20,6),ITBLE(50,2),DP(20),SCFS(20),C(20), 

&ZALPHA(20) ,IEND(6) ,DA(6) ,DIST(6) ,SEGN(6) ,DT(6) ,PEAK(6) ,ISG(6) , 

&NPU , NHD , NER , MAXNO , NCOMM , ICC , NCODE , TIME , KCODE , ICODE 

DATA  ZALPHA/IH1,1H2,1H3,1H4,IH5,1H6,1H7,1H8,1H9,1H0,1H  , 

&1H*,1H.,1H,,1H-,1H  ,IH  ,1H  , IH  , IH  / 

DATA  NCOMM/ 17/ 

DATA  CTBLE/ IHS , IHS , IHR, IHC , IHP , IHP, IHP, IHA, IHS , IHC , IHS , IHC , IHR, 
&1HR,1HE,IHS,1HF,33*1H  , 

&1HT,1HT,1HE,1H0,1HR,1HU,1HL,1HD,1HT,1H0,1HT,1H0,1H0,1H0,1HR, 

&1HE, 1HI,33*1H  , 

&2HAR , 2HOR , 2HCA , 2HMP , 2HIN , 2HNC , 2HOT , 2HD  , 2HOR, 2HMP , 2HOR , 2HMP , 

&2HUT , 2HUT , 2HRO , 2HDI , 2HNI , 33*2H  , 

&2HT  ,2HE  ,2HLL,2HUT,2HT  ,2HH  ,2H  H,2HHY,2HE  ,2HUT,2HE  ,2HUT, 

&2HE  ,2HE  ,2HR  , 2HME , 2HSH, 33*2H  , 


&2H 

,2HHY,2H  H,2HE  ,2HHY 

,2HHY,2HYD,2HD 

,2HRA,2HE  ,2HTR,2HE  , 

&2H 

,2HRE,2HAN,2HNT,2H 

,33*2H  , 

&2H 

,2HD  ,2HYD,2HHY,2HD 

, 2HD  , 2H  , 2H 

, 2HTI , 2HRA , 2HAV , 2HTR, 

&2H 

,2HSE,2HAL,2H  Y,2H 

,33*2H  , 

&2H 

,2H  ,2H  ,2HD  , 2H 

,2H  ,2H  ,2H 

,2HNG,2HTI,2HEL,2HAV, 

&2H 

,2HRV,2HYS,2HIE,2H 

,33*2H  , 

&8*2H  ,2H  C,2HNG,2H  T,2HEL,2H  , 2HOI , 2HIS , 2HLD, 34*2H  , 

&8*2H  ,2HUR,2H  C,2HIM,2H  T,2H  ,2HR  ,36*2H  , 

&8*2H  ,2HVE,2HUR,2HE  ,2HIM,38*2H  , 

&9*2H  ,2HVE,2H  ,2HE  ,38*2H'  / 


DATA  ITBLE/ 1,2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1 7, 33*1H 
&2, 310, 310, 310, 2, 1,2, 4, 100, 310, 100, 5, 4, 100, 2, 5, 0,33*1H  / 

END 


/ 


Data  files  for  MILHY2 


* 


*  NORTH  CREEK  -  TEXAS 

*  STORM  6 

* 

* 

START  RAINFALL  BEGINS  AT  15.25  HRS 

COMPUTE  HYD  ID=1  HYDNO=354  DT-.25  DA=23.7  CN=87  HT=355  L=1 1.745 

CUMULATIVE  RAINFALL  =  .0  .0  .01  .03  .13  .16  .27  .94 
1.11  1.15  1.19  1.2  1.22  1.46  1.48  1.5  1.53  1.54  1.55 
1.56  1.57  1.6  1.63  1.66  1.69  1.69  1.7  1.7  1.71  1.72 
1.73  1.74  1.75  1.76  1.77  1.78  1.78 

* 

* 


STORE  HYD 


PRINT  HYD 
PRINT  HYD 
PLOT  HYD 

* 


ID=2  HYD  NO=454  DT=».25  DA-23.7  FLOW  RATES- 
.0  .0  .0  .0  50.  15.  69.  122.  333.  417.  549.  680.  710. 
740,  825.  910.  1030.  1150.  1250.  1350.  1393.  1435. 
1478.  1520.  1490.  1460.  1430.  1400.  1320.  1240.  1160. 
1080.  1005.  930.  855.  780.  730.  680.  630.  580.  546. 


512. 

.  478 

• 

444. 

419. 

394. 

,  369. 

344. 

324. 

304. 

.  284 

• 

264. 

248. 

.  231 

• 

215. 

198. 

189. 

.  179. 

170. 

161. 

152, 

.  142 

• 

137. 

132. 

.  127 

• 

122. 

118. 

114. 

.  110. 

106. 

102. 

98. 

96. 

92 

.  89. 

86. 

83. 

80 

.  78. 

76. 

74. 

72.  70 

.  68. 

66. 

65. 

63. 

61 

.  59. 

57. 

56. 

55 

.  54. 

53. 

52. 

51.  50 

.  48. 

46. 

44. 

43. 

42 

.  41. 

40. 

39. 

38 

.  37. 

36. 

35. 

35.  35 

.  34. 

34. 

34. 

34. 

33 

.  33. 

33. 

33. 

32 

.  32. 

32. 

32. 
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15.25  15.25  0.25  9.0 
1 

60.  .25 
3 

23 

10  4  2 

.15  .15  .15  .15  .15  .15  .15  .15  .15  .15 

0.  .0000  .000 

.4300  .000  .4500  .000  .3600  .000 
6.940-7  .000  1.50-7  .000  4.44D-7  .000 

.4300  .4300  .4300  .4300  .4400  .4400  .3500  .3500  .3500  .3500 
.000 
6 

.242700  .27200  .30000  .33500  .37100  .43400 
-20.  -10.  -6.  -3.  -2.  -.4 
.000 

.29800  .32600  .35600  .38500  .41600  .45000 
-20.  -10.  -6.  -3.  -2.  -.8 
.000 

.25900  .27600  .29400  .31000  .33000  .3600 
-20.  -10.  -6.  -3.  -2.  -.35 
0.000 
67 

6  1  3 

.15  .15  .15  .15  .15  .15 

0.  .0000  .000 

.42500  .000  .4800  .000  .4800  .000 
7.20-6  .000  1.670-8  .000  1.670-8  .000 
.4200  .4600  .4600  .4600  .4600  .4600 
0.000 

7 

.017100  .10900  .12700  .14800  .1700  .21000  .42500 
-150.  -20.  -10.  -6.  -3.  -2.  -.3 
.000 

.31800  .37800  .40400  .42900  .45400  .47800  .4800 

-150.  -20.  -10.  -6.  -3.  -2.  -.3 

.000 

.31800  .37800  .40400  .42900  .45400  .47800  .4800 
-150.  -20.  -10.  -6.  -3.  -2.  -.3 
0.000 
10 

9  1  4 

.15  .15  .15  .15  .15  .15  .15  .15  .15 

0.0  0.000  0.000 

.3600  .000  .4800  .000  .4800  .000 
6.390-7  .000  1.670-7  .000  1.670-7  .000 
.3500  .4700  .4700  .4700  .4700  .4700  .4700  .4700  .4700 
0.000 
6 

.26100  .28600  .3100  .33600  .3500  .3600 
-20.  -10.  -6.  -3.  -2.  -.6 
.000 

.37800  .40400  .42900  .45400  .47800  .4800 
-20.  -10.  -6.  -3.  -2.  -.86 
.000 

.37800  .40400  .42900  .45400  .47800  .4800 

-20.  -10.  -6.  -3.  -2.  -.86 

0.000 


*  NORTH  CREEK  -  TEXAS 
\c 

*  STORM  6 
\c 

* 

\c 

* 

\c 

* 

\c 

START  RAINFALL  BEGINS  AT  15.25  HRS 

\c 

COMPUTE  HYD  ID=1  HYDNO-354  DT=.25  DA-23.7  CN-87  HT-355  L-1 1.745 

\c 

CUMULATIVE  RAINFALL  -  .0  .0  .01  .03  .13  .16  .27  .94 
\c 

1.11  1.15  1.19  1.2  1.22  1.46  1.48  1.5  1.53  1.54  1.55 
\c 

1.56  1.57  1.6  1.63  1.66  1.69  1.69  1.7  1.7  1.71  1.72 
\c 

1.73  1.74  1.75  1.76  1.77  1.78  1.78 
\c 

Shape  constant,  N  =  3.319 

Unit  peak  =  1567.3  cms 

INCREMENTAL  RUNOFF-Parameter  variability  included 
SD  of  detcap  0.000 

SD  of  saturated  soil  contentO.OOO  layer  1 

0,000  layer  2 
0,000  layer  3 

SD  of  suction  moisture  curveO.OOO  layer  1 

0.000  layer  2 
0.000  layer  3 

SD  of  sat  conductivityO.OOO  layer  1 

0,000  layer  2 
0.000  layer  3 

SD  of  Initial  water  contentO.OOO 

GENERATED  K-MOISTURE  CURVE 
Millington-Qulrk  Method 

Layer  1  Layer  2  Layer  3 

Moisture  Suction  Unsat  K  Moisture  Suction  Unsat  K  Moisture 
Sucti 

\con  Unsat  K 

0.243  -20.000  0.000000000033  0.298  -20.000  0.000000000014  0.259 


\cOO  0.000000000023 

0.253  -16.564  0.000000000146 

0.306 

-17.143 

0.000000000062 

-16.8 

\c73  0.000000000098 

0.263  -13.127  0.000000000373 

0.314 

-14.286 

0.000000000154 

-13.7 

\c46  0.000000000248 

0.273  -9.871  0.000000000789 

0.322 

-11.429 

0.000000000311 

-10.6 

\cl9  0.000000000515 

0.283  -8.432  0.000000001479 

0.330 

-9.467 

0.000000000562 

-9.0 

\c53  0.000000000953 

0.293  -6.994  0.000000002541 

0.338 

-8.400 

0.000000000937 

-7.8 

\c71  0.000000001613 

0.303  -5.733  0.000000004132 

0.346 

-7.333 

0.000000001466 

-6.6 

\c90  0.000000002568 

0.313  -4.870  0.000000006446 

0.354 

-6.267 

0.000000002195 

-5.5 

\c86  0.000000003931 

0.323  -4.007  0.000000009764 

0.362 

-5.379 

0.000000003191 

-4.5 

\c89  0.000000005880 

0.333  -3.144  0.000000014612 

0.370 

-4.552 

0.000000004545 

-3.5 

\c92  0.000000008738 

0.343  -2.767  0.000000021541 

0.378 

-3.724 

0.000000006409 

-2.8 

\c92  0.000000013010 

0.353  -2.487  0.000000030938 

0.386 

-2.968 

0.000000009050 

-2.6 

\c26  0.000000019126 

0.364  -2.208  0.000000043250 

0.394 

• 

-2.710 

0.000000012702 

-2.3 

\c61  0.000000027424 

0.374  -1.934  0.000000059123 

0.402 

-2.452 

0.000000017521 

-2.0 

\c95  0.000000038379 

0.384  -1.679  0.000000079470 

0.410 

-2.194 

0.000000023725 

-1.8 

\cl2  0.000000052724 

0.394  -1.423  0.000000105690 

0.418 

-1.929 

0.000000031631 

-1.5 

\cl9  0.000000071686 

0.404  -1.167  0.000000140181 

0.426 

-1.647 

0.000000041757 

-1.2 

\c27  0.000000097445 

0.414  -0.911  0.000000187554 

0.434 

-1.365 

0.000000054986 

\c35  0.000000134366 


0.424  -0.656  0.000000258405 

0.442 

-1.082  0.000000072941 

0.355 

-0.6 

\c42  0.000000193370 

0.434  -0.40  0.000000387204 

0.450 

-0.800  0.000000099113 

0.360 

-0. 

\c35  0.000000318548 

START  CONDITIONS 

Simulation  start  timel5.3hrs 

Precipitation  begins  at  15.3  and  ends  at  24.3 

Rainfall  data  time  increment  -  0.2500  hrs 

Time  increment  for  iteration  period  =  60.0  secs 

Maximum  evaporation  during  the  day  »  0.00000000  ms-1 
Surface  detention  capacity  =  0.0000  m 


INITIAL  SOIL  COLUMN  CONDITIONS 


SAT 

SAT  HYD 

CELL  DEPTH 

INITAL 

REL 

THETA 

COND 

NO 

THETA 

SAT 

m3/m3 

ms-1 

m 

m3/m3 

Layer 

1 

0.4350 

0.000000694000 

1 

0.0750 

0.4300 

0.989 

2 

0.2250 

0.4300 

0.989 

3 

0.3750 

0.4300 

0.989 

4 

0.5250 

0.4300 

0.989 

Layer 

2 

0.4510 

0.000000150000 

5 

0.6750 

0.4400 

0.976 

6 

0.8250 

0.4400 

0.976 

Layer 

3 

0.3610 

0.000000444000 

7 

0.9750 

0.3500 

0.970 

8 

1.1250 

0.3500 

0.970 

. 

9 

1.2750 

0.3500 

0.970 

10 

1.4250 

0.3500 

0.970 

SOIL  COLUMN  CONDITIONS  0.250  HRS  SINCE  RAIN  BEGAN 


Cell 

Depth  SWP 

Theta 

Hyd  cond 

Net  flux 

Rel  sat 

1 

0.0750  -0.5645 

0.4273 

0.000000304 

-0.00003174 

0.982 

2 

0.2250  -0.7034 

0.4216 

0.000000245 

-0.00006666 

0.969 

3 

0.3750  -2.2778 

0.3610 

0.000000040 

0.000095139 

0.830 

4 

0.5250  -2.3384 

0.3587 

0.000000037 

-0.00002351 

0.825 

5 

0.6750  -5.4902 

0.3610 

0.000000003 

0.000026605 

0.800 

6 

0.8250  -5.4902 

0.3610 

0.000000003 

0.000072302 

0.800 

7 

0.9750  -1.3438 

0.3415 

0.000000087 

-0.00006142 

0.946 

8 

1.1250  -0.9503 

0.3490 

0.000000132 

-0.00001637 

0.967 

9 

1.2750  -0.9044 

0.3499 

0.000000140 

-0.00000253 

0.969 

10 

1.4250  -0.9003 

0.3500 

0.000000141 

-0.00000025 

0.970 

-141- 


Balance  check  on  soil  column  water  status  »  -0.047726 

Balance  check  as  column  water  vol.  *  -8.6438310  % 

Cumulative  evaporation  =  0.00000000 

Cumulative  precipitation  =  0.0000 

Cumulative  infiltration  ”  0.000000 

Cumulative  drainage  =  0.000127 

Detention  capacity  exceeded 

Runoff  total  in  the  last  period  0.0000000  m 

Runoff  total  in  the  last  period  0.0000000  ins  0.250 


SOIL  COLUMN  CONDITIONS  0.500  HRS  SINCE  RAIN  BEGAN 


Cell 

Depth  SWP 

Theta 

Hyd  cond 

Net  flux 

Rel  sat 

1 

0.0750  -0.6125 

0.4255 

0.000000280 

-0.00001982 

0.978 

2 

0.2250  -0.8394 

0.4164 

0.000000207 

-0.00004191 

0.957 

3 

0.3750  -2.2778 

0.3610 

0.000000040 

0.000074599 

0.830 

4 

0.5250  -2.3996 

0.3565 

0.000000035 

-0.00002046 

0.819 

5 

0.6750  -5.4902 

0.3610 

0.000000003 

0.000024360 

0.800 

6 

0.8250  -5.4902 

0.3610 

0.000000003 

0.000051867 

0.800 

7 

0.9750  -1.6053 

0.3369 

0.000000066 

-0.00003634 

0.933 

8 

1.1250  -1.0445 

0.3473 

0.000000121 

-0.00001698 

0.962 

9 

1.2750  -0.9266 

0.3495 

0.000000136 

-0.00000541 

0.968 

10 

1.4250  -0.9043 

0.3499 

0.000000141 

-0.00000136 

0.969 

Balance  check  on  soil 

column  water  status  = 

-0.050252 

Balance  check  as  column  water  vol.  =*  -9.1409788  % 

Cumulative  evaporation  =  0.00000000 

Cumulative  precipitation  »  0.0003 

Cumulative  infiltration  =  0.000254 

Cumulative  drainage  =  0.000254 

Detention  capacity  exceeded 

Runoff  total  in  the  last  period  0.0000000  m 

Runoff  total  in  the  last  period  0.0000000  ins  0.500 


SOIL  COLUMN  CONDITIONS  0.750  HRS  SINCE  RAIN  BEGAN 


Cell 

Depth  SWP 

Theta 

Hyd  cond 

Net  flux 

Rel  sat 

1 

0.0750  -0.6282 

0.4250 

0.000000272 

-0.00000704 

0.977 

2 

0.2250  -0.9256 

0.4131 

0.000000185 

-0.00002670 

0.950 

3 

0.3750  -2.2778 

0.3610 

0.000000040 

0.000062898 

0.830 

4 

0.5250  -2.4531 

0.3546 

0.000000032 

-0.00001791 

0.815 

5 

0.6750  -5.4902 

0.3610 

0.000000003 

0.000022456 

0.800 

6 

0.8250  -5.4902 

0.3610 

0.000000003 

0.000041719 

0.800 

7 

0.9750  -1.7741 

0.3339 

0.000000055 

-0.00002538 

0.925 

8 

1.1250  -1.1332 

0.3457 

0.000000109 

-0.00001504 

0.958 

9  1.2750  -0.9603  0.3489  0.000000131  -0.00000675  0.966 

10  1.4250  -0.9152  0.3497  0.000000138  -0.00000265  0.969 

Balance  check  on  soli  column  water  status  =  -0.052308 

Balance  check  as  column  water  vol.  =  -9.5440872  % 

Cumulative  evaporation  =  0.00000000 

Cumulative  precipitation  =  0.0008 

Cumulative  infiltration  =  0.000762 

Cumulative  drainage  =  0.000380 

Detention  capacity  exceeded 

Runoff  total  in  the  last  period  0.0000000  m 

Runoff  total  in  the  last  period  0.0000000  ins  0.750 


SOIL  COLUMN  CONDITIONS  1.000  HRS  SINCE  RAIN  BEGAN 


Cell 

Depth  SWP 

Theta 

Hyd  cond 

Net  flux 

Rel  sat 

1 

0.0750  0.0000 

0.4344 

0.000000385 

-0.00009129 

0.999 

2 

0.2250  -0.9405 

0.4131 

0.000000182 

0.000057525 

0.950 

3 

0.3750  -2.2778 

0.3610 

0.000000040 

0.000060876 

0.830 

4 

0,5250  -2.5000 

0,3529 

0.000000031 

-0.00001583 

0.811 

5 

0,6750  -5.4902 

0.3610 

0.000000003 

0.000020908 

0.800 

6 

0.8250  -5,4902 

0.3610 

0.000000003 

0.000035737 

0.800 

7 

0,9750  -1.8947 

0.3317 

0.000000049 

-0.00001974 

0.919 

8 

1,1250  -1,2106 

0.3443 

0.000000100 

-0.00001302 

0.954 

9 

1.2750  -0.9988 

0.3482 

0.000000126 

-0.00000716 

0.964 

10 

1.4250  -0.9326 

0.3494 

0.000000135 

-0.00000371 

0.968 

Balance 

!  check  on  soil 

column  water  status  = 

-0.054123 

Balance  check  as  column  water 

vol.  =  -9. 

8668475  % 

Cumulative  evaporation 

=  0. 

00000000 

Cumulative  precipitation  =  0. 

0033 

Cumulative  infiltration  =  0. 

003165 

Cumulative  drainage 

=  0, 

000503 

Detention  capacity  exceeded 

Runoff  total  in  the  last  period  0.0000000  m 

Runoff  total  in  the  last  period  0.0000000  ins  1.000 


SOIL  COLUMN  CONDITIONS  1.250  HRS  SINCE  RAIN  BEGAN 


Cell 

Depth 

SWP 

Theta 

Hyd  cond 

Net  flux 

Rel  sat 

1 

0.0750 

-0.430 

0.4327 

0.000000372 

-0.00001958 

0.995 

2 

0.2250  -0.9091 

0.4140 

0.000000188 

0.000001018 

0.952 

3 

0.3750  -2.2778 

0.3610 

0.000000040 

0.000063627 

0.830 

4 

0.5250  -2.5418 

0.3514 

0.000000029 

-0.00001420 

0.808 

5 

0.6750  -5.4902 

0.3610 

0.000000003 

0.000019756 

0.800 

6 

0.8250  -5.4902 

0.3610 

0.000000003 

0.000031517 

7 

0.9750  -1.9899 

0.3300 

0.000000044 

-0.00001598 

8 

1.1250  -1.2781 

0.3430 

0.000000093 

-0.00001151 

9 

1.2750  -1.0384 

0.3474 

0.000000121 

-0.00000722 

10 

1.4250  -0.9552 

0.3490 

0.000000132 

-0.00000452 

0.800 

0.914 

0.950 

0.962 

0.967 


Balance  check  on  soil  column  water  status  =  -0.055876 

Balance  check  as  column  water  vol.  =  -10.2044863  % 


Cumulative  evaporation 
Cumulative  precipitation 
Cumulative  infiltration 
Cumulative  drainage 


0.00000000 

0.0041 

0.004065 

0.000622 


Detention  capacity  exceeded 

Runoff  total  in  the  last  period  0.0000000  m 

Runoff  total  in  the  last  period  0.0000000  ins  1.250 


SOIL  COLUMN  CONDITIONS  1.500  HRS  SINCE  RAIN  BEGAN 


Cell 

Depth  SWP 

Theta 

Hyd  cond 

Net  flux 

Rel  sat 

1 

0.0750  -0.404 

0.4347 

0.000000385 

0.000135718 

0.999 

2 

0.2250  -0.8589 

0.4159 

0.000000202 

-0.00000496 

0.956 

3 

0.3750  -2.2778 

0.3610 

0.000000040 

0.000069876 

0.830 

4 

0.5250  -2.5794 

0.3501 

0.000000028 

-0.00001278 

0.805 

5 

0.6750  -5.4902 

0.3610 

0.000000003 

0.000018742 

0.800 

6 

0.8250  -5.4902 

0.3610 

0.000000003 

0.000028204 

0.800 

7 

0.9750  -2.0678 

0.3285 

0.000000040 

-0.00001325 

0.910 

8 

1.1250  -1.3382 

0.3420 

0.000000088 

-0.00001027 

0.947 

9 

1.2750  -1.0780 

0.3467 

0.000000116 

-0.00000712 

0.960 

10 

1.4250  -0.9817 

0.3485 

0.000000128 

-0.00000507 

0.965 

Balance 

:  check  on  soil 

column  water  status  = 

-0.057849 

Balance  check  as  column  water  vol.  -  -10.5678383  % 


Cumulative  evaporation 
Cumulative  precipitation 
Cumulative  infiltration 
Cumulative  drainage 


0.00000000 

0.0069 

0.005997 

0.000739 


Detention  capacity  exceeded 

Runoff  total  in  the  last  period  0.0003078  m 

Runoff  total  in  the  last  period  0.0121178  ins  1.500 


SOIL  COLUMN  CONDITIONS  1.750  HRS  SINCE  RAIN  BEGAN 

Cell  Depth  SWP  Theta  Hyd  cond  Net  flux 
1  0.0750  0.0000  0.4339  0.000000386  -0.00008395 


Rel  sat 
0.997 


2  0.2250 

3  0.3750 

4  0.5250 

5  0.6750 

6  0.8250 

7  0.9750 

8  1.1250 

9  1.2750 

10  1.4250 


-0.8237 

-2.2778 

-2.6132 

-5.4902 

-5.4902 

-2.1311 

-1.3920 

-1.1167 

-1.0105 


0.4176 

0.3610 

0.3488 

0.3610 

0.3610 

0.3273 

0.3410 

0.3460 

0.3480 


0.000000212 

0.000000040 

0.000000027 

0.000000003 

0.000000003 

0.000000037 

0.000000083 

0.000000111 

0.000000125 


0.000035491 

0.000074357 

-0.00001153 

0.000017844 

0.000025822 

-0.00001152 

-0.00000924 

-0.00000693 

-0.00000541 


Balance  check  on  soil  column  water  status  =  -0.060018 

Balance  check  as  column  water  vol.  =  -10.9756570  % 

Cumulative  evaporation  =  0.00000000 

Cumulative  precipitation  =  0.0239 

Cumulative  infiltration  =  0.007703 

Cumulative  drainage  =  0.000853 

Detention  capacity  exceeded 

Runoff  total  in  the  last  period  0.0143848  m 

Runoff  total  in  the  last  period  0.5663301  ins  1.750 


SOIL  COLUMN  CONDITIONS  2.000  HRS  SINCE  RAIN  BEGAN 


Balance  check  on  soil  column  water  status  =  -0.062239 

Balance  check  as  column  water  vol.  =  -11.3912754  % 

Cumulative  evaporation  =  0.00000000 

Cumulative  precipitation  =  0.0282 

Cumulative  infiltration  =  0.009580 

Cumulative  drainage  =  0.000964 

Detention  capacity  exceeded 

Runoff  total  in  the  last  period  0.0023587  m 

Runoff  total  in  the  last  period  0.0928615  ins  2.000 


SOIL  COLUMN  CONDITIONS  2.250  HRS  SINCE  RAIN  BEGAN 


0.960 

0.830 

0.802 

0.800 

0.800 

0.907 

0.945 

0.958 

0.964 


Cell 

Depth 

SWP 

Theta 

Hyd  cond 

Net  flux 

Rel  sat 

1 

0.0750 

0.0000 

0.4344 

0.000000386 

-0.00008278 

0.999 

2 

0.2250 

-0.8054 

0.4183 

0.000000217 

0.000031749 

0.961 

3 

0.3750 

-2.2778 

0.3610 

0.000000040 

0.000076631 

0.830 

4 

0.5250 

-2.6438 

0.3478 

0.000000026 

-0.00001043 

0.799 

5 

0.6750 

-5.4902 

0.3610 

0.000000003 

0.000017046 

0.800 

6 

0.8250 

-5.4902 

0.3610 

0.000000003 

0.000023964 

0.800 

7 

0.9750 

-2.1858 

0.3262 

0.000000035 

-0.00001029 

0.904 

8 

1.1250 

-1.4406 

0.3401 

0.000000079 

-0.00000842 

0.942 

9 

1.2750 

-1.1543 

0.3453 

0.000000107 

-0.00000671 

0.957 

10 

1.4250 

-1.0408 

0,3474 

0.000000121 

-0.00000559 

0.962 

Cell 

Depth  SWP 

Theta 

Hyd  cond 

Net  flux 

Rel  sat 

1 

0.0750  -0.402 

0.4349 

0.000000386 

-0.00000096 

1.000 

2 

0.2250  -0.8285 

0.4171 

0.000000211 

-0.00001148 

0.959 

3 

0.3750  -2.2778 

0.3610 

0.000000040 

0.000073128 

0.830 

4 

0.5250  -2.6715 

0.3468 

0.000000025 

-0.00000945 

0.797 

5 

0.6750  -5.4902 

0.3610 

0.000000003 

0.000016334 

0.800 

6 

0.8250  -5.4902 

0.3610 

0.000000003 

0.000022338 

0.800 

7 

0.9750  -2.2348 

0.3252 

0.000000033 

-0.00000928 

0.901 

8 

1.1250  -1.4852 

0.3393 

0.000000075 

-0.00000775 

0.940 

9 

1.2750  -1.1906 

0.3447 

0.000000102 

-0.00000647 

0.955 

10 

1.4250  -1.0718 

0.3468 

0.000000117 

-0.00000565 

0.961 

Balance 

:  check  on  soil 

column  water  status  = 

-0.063968 

Balance 

1  check  as  column  water 

vol.  =  -11.7231504  Z 

Cumulative  evaporation 

=  0. 

00000000 

Cumulative  precipitation  =  0. 

0292 

Cumulative  infiltration  =  0. 

010700 

Cumulative  drainage 

=  0. 

001071 

Detention  capacity  exceeded 

Runoff 

total  in  the  last  period 

0.0000000  m 

Runoff 

total  in  the  last  period 

0.0000000  ins 

;  2.250 

SOIL  COLUMN  CONDITIONS 

2.500 

HRS  SINCE  RAIN 

BEGAN 

Cell 

Depth  SWP 

Theta 

Hyd  cond 

Net  flux 

Rel  sat 

1 

0.0750  -0.402 

0.4347 

0.000000386 

-0.00000290 

0.999 

2 

0.2250  -0.8510 

0.4162 

0.000000204 

-0.00000644 

0.957 

3 

0.3750  -2.2778 

0.3610 

0.000000040 

0.000069795 

0.830 

4 

0.5250  -2.6966 

0.3459 

0.000000024 

-0.00000859 

0.795 

5 

0.6750  -5.4902 

0.3610 

0.000000003 

0.000015697 

0.800 

6 

0.8250  -5.4902 

0.3610 

0.000000003 

0.000020902 

0.800 

7 

0.9750  -2.2792 

0.3244 

0.000000031 

-0.00000841 

0.899 

8 

1.1250  -1.5264 

0.3386 

0.000000071 

-0.00000720 

0.938 

9 

1.2750  -1.2256 

0.3440 

0.000000098 

-0.00000624 

0.953 

10 

1.4250  -1.1029 

0.3463 

0.000000113 

-0.00000563 

0.959 

Balance  check  on  soil 

column  water  status  = 

-0.065599 

Balance  check  as  column  water 

vol.  =  -12.0379328  % 

Cumulative  evaporation 

=  0. 

00000000 

Cumulative  precipitation  *  0. 

0302 

Cumulative  infiltration  =  0. 

011716 

Cumulative  drainage 

*  0. 

001174 

Detention  capacity  exceeded 

Runoff 

total  in  the  last  period 

0.0000000  m 

Runoff 

total  in  the  last  period 

0.0000000  ins 

;  2.500 
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SOIL  COLUMN  CONDITIONS 

2.750 

HRS  SINCE  RAIN 

BEGAN 

Cell 

Depth  SWP 

Theta 

Hyd  cond 

Net  flux 

Rel  sat 

1 

0.0750  -0.497 

0.4299 

0.000000338 

-0.00003988 

0.988 

2 

0.2250  -0.8778 

0.4151 

0.000000197 

-0.00001666 

0.954 

3 

0.3750  -2.2778 

0.3610 

0.000000040 

0.000065988 

0.830 

4 

0.5250  -2.7195 

0.3450 

0.000000023 

-0.00000781 

0.793 

5 

0.6750  -5.4902 

0.3610 

0.000000003 

0.000015126 

0.800 

6 

0.8250  -5.4902 

0.3610 

0.000000003 

0.000019630 

0.800 

7 

0.9750  -2.3194 

0.3236 

0.000000029 

-0.00000761 

0.896 

8 

1.1250  -1.5649 

0.3379 

0.000000069 

-0.00000675 

0.936 

9 

1.2750  -1.2595 

0.3434 

0.000000095 

-0.00000606 

0.951 

10 

1.4250  -1.1338 

0.3457 

0.000000109 

-0.00000556 

0.958 

Balance 

1  check  on  soil  ^ 

column  water  status  » 

-0.067155 

Balance  check  as  column  water 

Cumulative  evaporation  =  0. 

Cumulative  precipitation  =  0. 

Cumulative  infiltration  =  0. 

Cumulative  drainage  »  0. 

Detention  capacity  exceeded 
Runoff  total  in  the  last  period 

vol.  -  -12.3551773  % 

00000000 

0305 

011970 

001274 

0.0000000  m 

Runoff 

total  in  the  last  period 

0.0000000  ins 

!  2.750 

SOIL  COLUMN  CONDITIONS 

3.000 

HRS  SINCE  RAIN 

BEGAN 

Cell 

Depth  SWP 

Theta 

Hyd  cond 

Net  flux 

Rel  sat 

1 

0.0750  -0.5510 

0.4279 

0.000000311 

-0.00001778 

0.984 

2 

0.2250  -0.9206 

0.4134 

0.000000186 

-0.00001647 

0.950 

3 

0.3750  -2.2778 

0.3610 

0.000000040 

0.000060459 

0.830 

4 

0.5250  -2.7403 

0.3443 

0.000000022 

-0.00000712 

0.792 

5 

0.6750  -5.4902 

0.3610 

0.000000003 

0.000014611 

0.800 

6 

0.8250  -5.4902 

0.3610 

0.000000003 

0.000018499 

0.800 

7 

0.9750  -2.3558 

0.3228 

0.000000028 

-0.00000694 

0.894 

8 

1.1250  -1.6010 

0.3372 

0.000000066 

-0.00000635 

0.934 

9 

1.2750  -1.2924 

0.3428 

0.000000092 

-0.00000587 

0.950 

10 

1.4250  -1.1641 

0.3452 

0.000000105 

-0.00000546 

0.956 
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Balance  check  on  soil  column  water  status  =  -0.068606 

Balance  check  as  column  water  vol.  =  -12.6464275  X 

Cumulative  evaporation  =  0.00000000 

Cumulative  precipitation  =  0.0310 
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